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Introduction 


This  report  is  a FORTRAN  program  to  compute  the  Green's  functions  of  an 
infinite  plate.  The  Green's  function,  G. . (£,x,t),  is  defined  as  the  ith  com- 
ponent of  the  displacement  at  x due  to  anoint  force  of  step-function  time 
dependency  acting  at  in  the  jth  direction  initiated  at  t=0.  The  Green's 
function  is  the  fundamental  solution  of  the  transient  elastic  wave  propagation 
problem.  In  general,  the  displacement  field  u(^,x,t)  at  x due  to  a point  force 
of  arbitrary  time  dependence  acting  at  can  be  computed  by  a convolution 
integration;  i . e. , 


u.  U,x,t ) 


CO 

G.  . (j(,x  ,t)f  . ( T-t  )dt. 


(1 ) 


Here,  G. . is  the  time  derivative  of  G. . and  f.(t)  is  the  point  force  component 
of  arbi^ary  time  dependence  acting  inJthe  jtA  direction  (summation  over 
repeated  indices  is  used).  Displacement  due  to  point  dipoles  or  couple  forces 
can  be  represented  by  the  spatial  derivatives  of  G. ..  Displacement  produced  by 
a dynamic  force  distributed  over  a finite  area  can1^lso  be  computed  by  numerical 
integration  using  the  Green's  function  as  the  kernel  of  the  integral  over  the 
finite  area. 


The  basic  formulation  of  the  problem  and  derivation  of  the  solution  for — 
mulas  were  reported  in  (1).  The  method  used  is  called  "ray  theory"  in  the 
seismological  literature.  Our  derivation  was  based  upon  John  Willis'  new 
Fourier  Inversion  method  (2).  Similar  computation  results  can  be  found  in 
Reference  (3),  (4),  and  (5). 

The  program  was  originally  written  for  the  analysis  of  acoustic  emission 
signals.  Acoustic  emission  is  a nondestructive  testing  and  monitoring  technique 
in  which  the  detection  of  the  transient  stress  wave  produced  by  localized  defor — 
mation  or  cracks  are  used  to  locate,  and  to  assess  the  criticality  of  the 
defects.  The  theoretical  computation  is  an  important  link  to  predict  how  the 
acoustic  emission  waveform  evolves  through  the  structure.  How  the  theoretical 
computations  are  used  in  the  study  of  acoustic  emission  can  be  found  in 
References  (1),  (6),  (7),  and  (8). 

This  computer  program  is  made  available  mainly  for  its  application  to 
calibrate  acoustic  emission  systems  and  sensors.  By  making  the  present  program 
available,  duplication  of  efforts  can  be  avoided,  errors  in  theory  can  be 
checked,  and  experimental  results  can  be  reproduced  and  verified. 

How  to  Use  the  Program 


I.  The  easiest  way  to  run  this  program  is  by  setting  up  nondimensional 
parameters  and  calling  the  subroutine  GREENFCT.  First,  the  x-y-z  coordinate.-' 


are  selected  by  choosing  an  x-y  plane  parallel  to  the  plate  and  the  origin  at 
the  center  of  the  plate  directly  underneath  the  source,  i.e.,  the  source  is 
always  located  at  £=(0,0,0. 5);  the  x axis  is  aligned  in  the  direction  from  the 
source  pointing  to  the  detector.  The  required  nondimensionalized  input 
parameters  for  the  subroutine  GREENFCT  are: 

1 . ALPHA  = shear  wave  speed/longitudinal  wave  speed; 

2.  XD  = x-coordinate  of  the  detector;  actual  length/thickness  of  the  plate; 

3.  ZD  = z-coordinate  of  the  detector;  actual  length/thickness  of  the  plate. 

ZD  = 0.5  if  the  source  and  the  detector  are  on  the  same  side  of  the 

plate ; 

ZD  = -0.5  if  they  are  on  opposite  sides  of  the  plate. 

(Note:  the  y-coordinate  of  the  detector  is  always  zero  because  of  the  way  the 

coordinate  system  is  chosen.) 

M.  INDEX  = ij ; subscript  of  Green's  function,  integer  number  11,  12,  13,  21, 
22,  23,  31,  32,  or  33-  (111,  112,  113,  etc.  for  force  dipoles.) 

5.  TDELTA  = sampling  time  interval  in  terms  of  nondimensionalized  time  unit 

= actual  sampling  time  interval  * shear  wave  speed/thickness  of  the 
plate;  and 

6.  NP0INT  = Total  number  of  sampling  points  to  be  computed;  must  be  an 

integer. 

The  subroutine  will  return  a double  precision  array  DISPL  of  dimension 
(NP0INT)  which  corresponds  to  the  desired  G (£,x,t)  or  G„  (£,x,t)  sampled  at 
equal  time  intervals.  However,  the  displacements  corresponding  to  dipoles  or 
couples  (i.e.,  G ^)  are  for  linear  ramp  time  dependency  input. 

Dif ferentitatingim?ie  returned  DISPL  with  respect  to  time  once  will  give  the 
proper  Green's  function  due  to  step  time  dependency  input. 

A simple  program  calls  GREENFCT  for  given  nondimensionalized  input  param- 
ters  together  with  the  output  results  is  included  in  Appendix  A;  a program  that 
prompts  the  user  for  input  parameters  of  arbitrary  physical  units  and  returns 
displacement  in  physical  units  is  included  in  Appendix  B.  The  complete  listing 
of  the  subroutines  is  in  Appexdix  C. 

II.  Another  way  to  use  this  program  is  by  calling  GREENSUB  which  is  basically  a 
simple  function:  for  given  one  time  T,  it  returns  a displacement  DISPL. 

However,  there  are  two  preparatory  steps  that  must  be  taken  in  the  main  program: 

(1)  RAYTIME  (NRAY , 3),  the  time  of  arrival  table,  TA(NRAY)  and  CN (NRAY , 3), 
two  working  arrays,  should  be  dimensioned. 

(2)  The  subroutines  INIT  and  TIMEARRI  must  be  called  first. 

The  subroutine  GREENFCT  may  serve  as  an  example  of  how  to  use  the  sub- 
routine GREENSUB. 

III.  In  addition  to  the  two  subroutines  mentioned  above,  there  are  three  sub- 
routines which  may  be  called  independently  and  the  user  may  find  them  useful  in 
checking  experimental  or  theoretical  results. 


(1)  EPIDIS  computes  the  vertical  displacement  at  the  epicenter  for  given 
ALPHA  and  T. 

(2)  SFWAVE  computes  the  surface  wave,  which  corresponds  to  the  Green's  func- 
tion of  a semi-infinite  space. 

(3)  TIMEARRI  computes  the  time  of  arrival  for  various  ray  paths. 

Comments  in  these  three  subroutines  serve  as  instructions  about  how  to  use 
these  subroutines. 

Some  Remarks 


I.  There  are  two  machine  dependent  constants  which  are  used  in  the  two  numeri- 
cal integration  routines  DGLQI  and  GLI5T.  EPMACH  is  the  largest  relative 
spacing.  UFLOW  is  the  smallest  positive  magnitude.  See  Reference  (9)  for 
details  about  how  to  set  these  constants  for  different  computers. 

II.  The  program  uses  double  precision  complex  arithmetic.  In  order  to  run  on 
those  FORTRAN  compilers  lacking  double  precision  complex  function  libraries,  the 
user  may  have  to  modify  the  program  to  single  precision  arithmetic. 

III.  The  ray  method  is  an  exact  solution  for  the  Green's  function  in  the  time 
domain  produced  by  summing  the  contribution  from  successive  arrivals  of 
reflected  rays.  The  numerical  computation  is  accurate  if  the  source  and  detec- 
tor are  near  each  other  (say  less  than  ten  plate  thicknesses)  and  if  the  maximum 
time  is  less  than  ten  dimensionless  units.  If  the  sensor  and  detector  are  far 
apart  or  the  maximum  time  is  large,  the  number  of  rays  arriving  at  approximately 
the  same  time  may  be  too  large  to  compute  in  a reasonable  time  and  the  accumu- 
lated error  may  grow. 

IV.  The  current  program  is  limited  to  the  test  configuration  where  the  sensor 
and  detector  are  on  the  surface  of  the  plate.  The  displacement  computed  for 
force  dipole  (i.e.,  INDEX  = 111,  112,  etc.)  are  for  point  dipole  of  linear  ramp 
time  dependency.  Differentiation  with  respect  to  time  once  of  the  displacement 
will  produce  proper  Green's  functions  due  to  step  time  dependency  (See  Reference 
10). 
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Appendix  A 


Sample  Program  and  Result 


^★★★★★★★★★★★★★★★★★★★★★★★★**************************************** 

C 


C 

c 

c 

c 

c 

c 


c 

c 


c 


PROGRAM  TSGREEN 

TEST  PROGRAM  TO  COMPUTE  GREEN'S  FUNCTION  OF  A PLATE 

DOUBLE  PRECISION  DISPLC<256) 

URITE(5,*) 'TYPE  INDEX, ALPHA, XD,  AND  ZD' 

READ(5,*)  INDEX, ALPHA, XD, ZD 

INDEX=33 

ALPHA=0. 59902 

XD=1 .0 

ZD=-0.5 

WRITE(5,*)  'TYPE  DELTAT  AND  NUMBER  OF  POINTS' 

READ(5,*)  TDELTA,NPOINT 

TDELTA=0.01 

NPOINT=256 


CALL  GREENFCTCALPHA, XD, ZD, INDEX, TDELTA,NPOINT,DISPLC) 

DO  1000  1=1,256 

TIME  = 0.01*REAL(I ) 

WRITEC6, 1001 ) TIME,DISPLC( I ) 

1000  CONTINUE 

1001  FORMAT(2X, 'TIME  = ',F8.3,', 

STOP 
END 


DISPLACEMENT  = \D21.14) 


TIME 

= 

0.010, 

DISPLACEMENT 

= 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.020, 

DISPLACEMENT 

= 

0.00000000000000D+00 

TIME 

= 

0.030, 

DISPLACEMENT 

= 

0.00000000000000D+00 

TIME 

= 

0.040, 

DISPLACEMENT 

= 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.050, 

DISPLACEMENT 

= 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.060, 

DISPLACEMENT 

= 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.070, 

DISPLACEMENT 

= 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.080, 

DISPLACEMENT 

= 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.090, 

DISPLACEMENT 

= 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.100, 

DISPLACEMENT 

= 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.110, 

DISPLACEMENT 

= 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.120, 

DISPLACEMENT 

= 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.130, 

DISPLACEMENT 

= 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.140, 

DISPLACEMENT 

= 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.150, 

DISPLACEMENT 

= 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.160, 

DISPLACEMENT 

= 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.170, 

DISPLACEMENT 

= 

O.OOOOOCOOOOOOOOD+OO 

TIME 

= 

0.180, 

DISPLACEMENT 

= 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.190, 

DISPLACEMENT 

= 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.200, 

DISPLACEMENT 

= 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.210, 

DISPLACEMENT 

= 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.220, 

DISPLACEMENT 

= 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.230, 

DISPLACEMENT 

= 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.240, 

DISPLACEMENT 

= 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.250, 

DISPLACEMENT 

= 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.260, 

DISPLACEMENT 

= 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.270, 

DISPLACEMENT 

= 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.280, 

DISPLACEMENT 

= 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.290, 

DISPLACEMENT 

= 

O.OOOOOOOOOOOOOOD+OO 

TIME 

- 

0.300, 

DISPLACEMENT 

= 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.310, 

DISPLACEMENT 

= 

0.000000000000000+00 

TIME 

= 

0.320, 

DISPLACEMENT 

r 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.330, 

DISPLACEMENT 

= 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.340, 

DISPLACEMENT 

= 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.350, 

DISPLACEMENT 

= 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.360, 

DISPLACEMENT 

= 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.370, 

DISPLACEMENT 

= 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.380, 

DISPLACEMENT 

= 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.390, 

DISPLACEMENT 

= 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.400, 

DISPLACEMENT 

= 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.410, 

DISPLACEMENT 

= 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.420, 

DISPLACEMENT 

= 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.430, 

DISPLACEMENT 

= 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.440, 

DISPLACEMENT 

= 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.450, 

DISPLACEMENT 

= 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.460, 

DISPLACEMENT 

= 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.470, 

DISPLACEMENT 

= 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.480, 

DISPLACEMENT 

= 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.490, 

DISPLACEMENT 

= 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.500, 

DISPLACEMENT 

= 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.510, 

DISPLACEMENT 

= 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.520, 

DISPLACEMENT 

= 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.530, 

DISPLACEMENT 

= 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.540, 

DISPLACEMENT 

= 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.550, 

DISPLACEMENT  = 

0.00000000000000D+00 

TIME 

= 

0.560, 

DISPLACEMENT  = 

0.00000000000000D+00 

TIME 

= 

0.570, 

DISPLACEMENT  = 

0.00000000000000D+00 

TIME 

= 

0.580, 

DISPLACEMENT  = 

0.00000000000000D+00 

TIME 

= 

0.590, 

DISPLACEMENT  = 

0.00000000000000D+00 

TIME 

= 

0.600, 

DISPLACEMENT  = 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.610, 

DISPLACEMENT  = 

0.00000000000000D+00 

TIME 

= 

0.620, 

DISPLACEMENT  = 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.630, 

DISPLACEMENT  = 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.640, 

DISPLACEMENT  = 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.650, 

DISPLACEMENT  = 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.660, 

DISPLACEMENT  = 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.670, 

DISPLACEMENT  = 

O.OOOOOOOOOOOOOOD+OO 

TIME 

- 

0.680, 

DISPLACEMENT  = 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.690, 

DISPLACEMENT  = 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.700, 

DISPLACEMENT  = 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.710, 

DISPLACEMENT  = 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.720, 

DISPLACEMENT  = 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.730, 

DISPLACEMENT  = 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.740, 

DISPLACEMENT  = 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.750, 

DISPLACEMENT  = 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.760, 

DISPLACEMENT  = 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.770, 

DISPLACEMENT  = 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.780, 

DISPLACEMENT  = 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.790, 

DISPLACEMENT  = 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.800, 

DISPLACEMENT  = 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.810, 

DISPLACEMENT  = 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.820, 

DISPLACEMENT  = 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.830, 

DISPLACEMENT  = 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.840, 

DISPLACEMENT  = 

O.OOOOOOOOOOOOOOD+OO 

TIME 

= 

0.850, 

DISPLACEMENT  = 

0.11 07506769871 2D+00 

TIME 

= 

0.860, 

DISPLACEMENT  = 

0.11 0986061 56782D+00 

TIME 

= 

0.870, 

DISPLACEMENT  = 

0.111 1 8024981 467D+00 

TIME 

= 

0.880, 

DISPLACEMENT  = 

0.11 133277774482D+00 

TIME 

= 

0.890, 

DISPLACEMENT  = 

0.11 14431 8482884D+00 

TIME 

= 

0.900, 

DISPLACEMENT  = 

0.1 1 151 1026521 21D+00 

TIME 

= 

0.910, 

DISPLACEMENT  = 

0.11 1535876944 15D+00 

TIME 

= 

0.920, 

DISPLACEMENT  = 

0.11 15173360582 1D+00 

TIME 

= 

0.930, 

DISPLACEMENT  = 

0.1 11455035416960+00 

TIME 

= 

0.940, 

DISPLACEMENT  = 

0.11 134864281 672D+00 

TIME 

= 

0.950, 

DISPLACEMENT  = 

0.1111 9787020707D+00 

TIME 

= 

0.960, 

DISPLACEMENT  = 

0.11 100247472 178D+00 

TIME 

= 

0.970, 

DISPLACEMENT  = 

0. 1 1076226552545D+00 

TIME 

= 

0.980, 

DISPLACEMENT  = 

0.1 104771 076421 9D+00 

TIME 

= 

0.990, 

DISPLACEMENT  = 

0. 1 1 014692394 960D+00 

TIME 

= 

1.000, 

DISPLACEMENT  = 

0.1 09771 71 429770D+00 

TIME 

= 

1.010, 

DISPLACEMENT  = 

0. 10935 1 509141 28D+00 

TIME 

= 

1.020, 

DISPLACEMENT  = 

0.1 08886450335 15D+00 

TIME 

= 

1.030, 

DISPLACEMENT  = 

0 . 1 0837678979646D+00 

TIME 

= 

1.040, 

DISPLACEMENT  = 

0 . 1 0782270845288D+00 

TIME 

= 

1.050, 

DISPLACEMENT  = 

0.1 0722461 062043D+00 

TIME 

= 

1.060, 

DISPLACEMENT  = 

0 . 1 0658293623297D+00 

TIME 

= 

1.070, 

DISPLACEMENT  = 

0 . 1 0589820949205D+00 

TIME 

= 

1.080, 

DISPLACEMENT  = 

0. 105171 1121 083 1 D+00 

TIME 

= 

1.090, 

DISPLACEMENT  = 

0 . 1 0440220346099D+00 

TIME 

= 

1.100, 

DISPLACEMENT  = 

0 . 1 0359233809357D+00 

TIME 

= 

1.110, 

DISPLACEMENT  = 

0 . 1 0274239332029D+00 

TIME 

- 

1.120, 

DISPLACEMENT  = 

0. 10185342388500D+00 

TIME 

= 

1.130, 

DISPLACEMENT  = 

0.1 009263 1946377D+00 

TIME 

= 

1.140, 

DISPLACEMENT  = 

0. 999623 13600975D- 01 

TIME 

= 

1.150, 

DISPLACEMENT  = 

0.9896264 7939981D- 01 

TIME 

= 

1.160, 

DISPLACEMENT  = 

0 . 9792875478595 1 D ■ 01 

TIME 

= 

1.170, 

DISPLACEMENT  = 

0. 96861 858478749D- 01 

TIME 

= 

1.180, 

DISPLACEMENT  = 

0. 957635660 16467D- 01 

TIME 

= 

1.190, 

DISPLACEMENT  = 

0 . 94635474207730D - 01 

TIME 

= 

1.200, 

DISPLACEMENT  = 

0 . 93479376354866D - 01 

TIME 

= 

1.210, 

DISPLACEMENT  = 

0 . 9229681 6327325D -01 

TIME 

= 

1.220, 

DISPLACEMENT  = 

0.91 089747924766D - 0 1 

TIME 

= 

1.230, 

DISPLACEMENT  = 

0.89860093569871D-01 

TIME 

= 

1.240, 

DISPLACEMENT  = 

0 .88609970379555D - 01 

TIME 

= 

1.250, 

DISPLACEMENT  = 

0. 87341 207961 754D- 01 

TIME 

= 

1.260, 

DISPLACEMENT  = 

0 . 86056060561 906D - 01 

TIME 

= 

1.270, 

DISPLACEMENT  = 

0. 84 7567281 27925D- 01 

TIME 

= 

1.280, 

DISPLACEMENT  = 

0. 834455961 72436D- 01 

TIME 

= 

1.290, 

DISPLACEMENT  = 

0 . 821 24730252424D - 01 

TIME 

= 

1.300, 

DISPLACEMENT  = 

0 . 80796622274542D - 01 

TIME 

= 

1.310, 

DISPLACEMENT  = 

0 . 79463686723 1 94D - 0 1 

TIME 

= 

1.320, 

DISPLACEMENT  = 

0.781 28380592089D - 01 

TIME 

= 

1.330, 

DISPLACEMENT  = 

0 . 76793325943265D - 01 

TIME 

= 

1.340, 

DISPLACEMENT  = 

0 . 75460795794622D - 01 

TIME 

= 

1.350, 

DISPLACEMENT  = 

0. 74133473 644589D- 01 

TIME 

= 

1.360, 

DISPLACEMENT  = 

0. 7281 393901 21 66D- 01 

TIME 

= 

1.370, 

DISPLACEMENT  = 

0. 71 50491 4223321 D- 01 

TIME 

= 

1.380, 

DISPLACEMENT  = 

0.70208761 123235D-01 

TIME 

= 

1.390, 

DISPLACEMENT  = 

0.68928225721454D-01 

TIME 

= 

1.400, 

DISPLACEMENT  = 

0 . 63777644085523D - 0 1 

TIME 

= 

1.410, 

DISPLACEMENT  = 

-0. 312054998264 77D-01 

TIME 

= 

1.420, 

DISPLACEMENT  = 

0. 367644440973 72D+00 

TIME 

= 

1.430, 

DISPLACEMENT  = 

0.391 058631 5 1275D+00 

TIME 

= 

1.440, 

DISPLACEMENT  = 

0. 40949206785671 D+00 

TIME 

= 

1.450, 

DISPLACEMENT  = 

0. 4258496636743 1D+00 

TIME 

= 

1.460, 

DISPLACEMENT  = 

0.441 02467848492D+00 

TIME 

= 

1.470, 

DISPLACEMENT  = 

0 . 45540362930489D+00 

TIME 

= 

1.480, 

DISPLACEMENT  = 

0. 4691 8701 3091 17D+00 

TIME 

= 

1.490, 

DISPLACEMENT  = 

0. 4824881 0354954D+00 

TIME 

= 

1.500, 

DISPLACEMENT  = 

0.4953792919241 1D+00 

TIME 

= 

1.510, 

DISPLACEMENT  = 

0 . 50790390548020D+00 

TIME 

= 

1.520, 

DISPLACEMENT  = 

0.52009143324477D+00 

TIME 

= 

1.530, 

DISPLACEMENT  = 

0.53 1961 91 904305D+00 

TIME 

= 

1.540, 

DISPLACEMENT  = 

0.543528255164190+00 

TIME 

= 

1.550, 

DISPLACEMENT  = 

0. 554802671341 78D+00 

TIME 

= 

1.560, 

DISPLACEMENT  = 

0. 565791 47374873D+00 

TIME 

= 

1.570, 

DISPLACEMENT  = 

0 . 57650038247907D+00 

TIME 

= 

1.580, 

DISPLACEMENT  = 

0. 5869330375 1993D+00 

TIME 

= 

1.590, 

DISPLACEMENT  = 

0. 597095366361 70D+00 

TIME 

= 

1.600, 

DISPLACEMENT  = 

0.6069898995 10390+00 

TIME 

= 

1.610, 

DISPLACEMENT  = 

0 .61661 997278485D+00 

TIME 

= 

1.620, 

DISPLACEMENT  = 

0.6259879 72 75653D+00 

TIME 

= 

1.630, 

DISPLACEMENT  = 

0 . 635 098976464060+00 

TIME 

= 

1.640, 

DISPLACEMENT  = 

0. 643955450564 72D+00 

TIME 

= 

1.650, 

DISPLACEMENT 

= 

0. 6525608866934 7D+00 

TIME 

= 

1.660, 

DISPLACEMENT 

= 

0.66091 81 3088975D+00 

TIME 

= 

1.670, 

DISPLACEMENT 

= 

0.6690325571 71 67D+00 

TIME 

= 

1.680, 

DISPLACEMENT 

= 

0. 67690728626371 D+00 

TIME 

= 

1.690, 

DISPLACEMENT 

= 

0. 684546391 70679D+00 

TIME 

= 

1.700, 

DISPLACEMENT 

= 

0 . 691 95338607382D+00 

TIME 

= 

1.710, 

DISPLACEMENT 

= 

0.6991 3401 239547D+00 

TIME 

= 

1.720, 

DISPLACEMENT 

= 

0.706091 992259960+00 

TIME 

= 

1.730, 

DISPLACEMENT 

= 

0.71 2831 858301 27D+00 

TIME 

= 

1.740, 

DISPLACEMENT 

= 

0.71 9358221 54785D+00 

TIME 

= 

1.750, 

DISPLACEMENT 

= 

0. 725675 16353530D+00 

TIME 

= 

1.760, 

DISPLACEMENT 

= 

0.731788611 32 189D+00 

TIME 

= 

1.770, 

DISPLACEMENT 

= 

0.73 7702708981 97D+00 

TIME 

= 

1.780, 

DISPLACEMENT 

= 

0 . 7434222332603 1 D+00 

TIME 

= 

1.790, 

DISPLACEMENT 

= 

0. 74895 144993679D+00 

TIME 

= 

1.800, 

DISPLACEMENT 

= 

0.754296 19306931 D+00 

TIME 

= 

1.810, 

DISPLACEMENT 

= 

0 . 75946069053969D+00 

TIME 

= 

1.820, 

DISPLACEMENT 

= 

0 . 764449680843430+00 

TIME 

= 

1.830, 

DISPLACEMENT 

= 

0 . 7692674 1601 862D+00 

TIME 

= 

1.840, 

DISPLACEMENT 

= 

0.77391947485854D+00 

TIME 

= 

1.850, 

DISPLACEMENT 

= 

0.7784099991 1627D+00 

TIME 

= 

1.860, 

DISPLACEMENT 

= 

0 . 78274354050484D+00 

TIME 

= 

1.870, 

DISPLACEMENT 

= 

0. 786924 19559247D+00 

TIME 

= 

1.880, 

DISPLACEMENT 

= 

0. 7909571 8088098D+00 

TIME 

= 

1.890, 

DISPLACEMENT 

= 

0 . 79484643653490D+00 

TIME 

= 

1.900, 

DISPLACEMENT 

= 

0.86806930229072D+00 

TIME 

= 

1.910, 

DISPLACEMENT 

= 

0 . 87369508356078D+00 

TIME 

= 

1.920, 

DISPLACEMENT 

= 

0 . 879233 1 4459302D+00 

TIME 

= 

1.930, 

DISPLACEMENT 

= 

0 . 884687498388460+00 

TIME 

= 

1.940, 

DISPLACEMENT 

= 

0.89006261 168853D+00 

TIME 

= 

1.950, 

DISPLACEMENT 

= 

0 . 89536286267026D+00 

TIME 

= 

1.960, 

DISPLACEMENT 

= 

0. 900592044445 14D+00 

TIME 

= 

1.970, 

DISPLACEMENT 

= 

0.905755352541 14D+00 

TIME 

= 

1.980, 

DISPLACEMENT 

= 

0.91 085639280247D+00 

TIME 

= 

1.990, 

DISPLACEMENT 

= 

0.91 58991 8253638D+00 

TIME 

= 

2.000, 

DISPLACEMENT 

= 

0.9208871 7546076D+00 

TIME 

= 

2.010, 

DISPLACEMENT 

= 

0 . 925825 1 5984505D+00 

TIME 

= 

2.020, 

DISPLACEMENT 

= 

0.93071 640538221 D+00 

TIME 

= 

2.030, 

DISPLACEMENT 

= 

0.935564571 71 401 D+00 

TIME 

= 

2.040, 

DISPLACEMENT 

= 

0. 9403727744075 7D+00 

TIME 

= 

2.050, 

DISPLACEMENT 

= 

0.9451 454 1664089D+00 

TIME 

= 

2.060, 

DISPLACEMENT 

= 

0 . 94988543954263D+00 

TIME 

= 

2.070, 

DISPLACEMENT 

= 

0.9545961 61 31 288D+00 

TIME 

= 

2.080, 

DISPLACEMENT 

= 

0.95928037273261 D+00 

TIME 

= 

2.090, 

DISPLACEMENT 

= 

0.9639421 23260 14D+00 

TIME 

= 

2.100, 

DISPLACEMENT 

= 

0.96858404264596D+00 

TIME 

= 

2.110, 

DISPLACEMENT 

= 

0.97320913182266D+00 

TIME 

= 

2.120, 

DISPLACEMENT 

= 

0. 97781 987758092D+00 

TIME 

= 

2.130, 

DISPLACEMENT 

= 

0. 98242001 202845D+00 

TIME 

= 

2.140, 

DISPLACEMENT 

= 

0.98701 187638821 D+00 

TIME 

= 

2.150, 

DISPLACEMENT 

= 

0.991 5981 8288714D+00 

TIME 

= 

2.160, 

DISPLACEMENT 

= 

0 . 9961 81 57630830D+00 

TIME 

= 

2.170, 

DISPLACEMENT 

= 

0.1 0007641 9879 12D+01 

TIME 

= 

2.180, 

DISPLACEMENT 

= 

0. 10053494399595D+01 

TIME 

= 

2.190, 

DISPLACEMENT 

= 

0.1 0099393 168333D+01 

TIME 

= 

2.200, 

DISPLACEMENT  = 

0. 10145362237750D+01 

TIME 

= 

2.210, 

DISPLACEMENT  = 

0.10191420575858D+01 

TIME 

= 

2.220, 

DISPLACEMENT  = 

0.1 02375 99774067D+01 

TIME 

= 

2.230, 

DISPLACEMENT  = 

0.1 028391 771 2906D+01 

TIME 

= 

2.240, 

DISPLACEMENT  = 

0.1 03303961 40693D+01 

TIME 

= 

2.250, 

DISPLACEMENT  = 

0.1 037705 1845379D+01 

TIME 

= 

2.260, 

DISPLACEMENT  = 

0.1 042391 4489334D+01 

TIME 

= 

2.270, 

DISPLACEMENT  = 

0.1 047099993231 4D+01 

TIME 

= 

2.280, 

DISPLACEMENT  = 

0. 10518328029171D+01 

TIME 

= 

2.290, 

DISPLACEMENT  = 

0. 10565913655695D+01 

TIME 

= 

2.300, 

DISPLACEMENT  = 

0.106137849050690+01 

TIME 

= 

2.310, 

DISPLACEMENT  = 

0.1 1476631 005725D+01 

TIME 

= 

2.320, 

DISPLACEMENT  = 

0. 1 1533185695767D+01 

TIME 

= 

2.330, 

DISPLACEMENT  = 

0. 1 1590500263745D+01 

TIME 

= 

2.340, 

DISPLACEMENT  = 

0.11648631801910D+01 

TIME 

= 

2.350, 

DISPLACEMENT  = 

0.11 707621 604698D+01 

TIME 

= 

2.360, 

DISPLACEMENT  = 

0.11 7675 16974758D+01 

TIME 

= 

2.370, 

DISPLACEMENT  = 

0.1 182836581 1738D+01 

TIME 

= 

2.380, 

DISPLACEMENT  = 

0. 11890210681560D+01 

TIME 

= 

2.390, 

DISPLACEMENT  = 

0.1 19531 12520959D+01 

TIME 

= 

2.400, 

DISPLACEMENT  = 

0. 120171 15269290D+01 

TIME 

= 

2.410, 

DISPLACEMENT  = 

0. 12082269375580D+01 

TIME 

= 

2.420, 

DISPLACEMENT  = 

0. 12148619575064D+01 

TIME 

= 

2.430, 

DISPLACEMENT  = 

0.1 221 62303280 19D+01 

TIME 

= 

2.440, 

DISPLACEMENT  = 

0.1 2285 147870224D+01 

TIME 

= 

2.450, 

DISPLACEMENT  = 

0. 12355425415795D+01 

TIME 

= 

2.460, 

DISPLACEMENT  = 

0.124271099911260+01 

TIME 

= 

2.470, 

DISPLACEMENT  = 

0 . 1 2500269925297D+01 

TIME 

= 

2.480, 

DISPLACEMENT  = 

0.1 2574953844 797D+01 

TIME 

= 

2.490, 

DISPLACEMENT  = 

0.12651217887352D+01 

TIME 

= 

2.500, 

DISPLACEMENT  = 

0. 127291 11430751D+01 

TIME 

= 

2.510, 

DISPLACEMENT  = 

0. 12808706969625D+01 

TIME 

= 

2.520, 

DISPLACEMENT  = 

0. 12890055553986D+01 

TIME 

= 

2.530, 

DISPLACEMENT  = 

0.1 29732 16346622D+01 

TIME 

= 

2.540, 

DISPLACEMENT  = 

0. 1305824 1075455D+01 

TIME 

= 

2.550, 

DISPLACEMENT  = 

0. 13145206649727D+01 

TIME 

= 

2.560, 

DISPLACEMENT  = 

0.1 32341 665 14682D+01 

TIME 

= 

2.080, 

DISPLACEMENT  = 

0.959280372 73261 D+00 

TIME 

= 

2.090, 

DISPLACEMENT  = 

0 . 9639421 232601 4D+00 

TIME 

= 

2.100, 

DISPLACEMENT  = 

0 . 96858404264596D+00 

TIME 

= 

2.110, 

DISPLACEMENT  = 

0.973209131822660+00 

TIME 

= 

2.120, 

DISPLACEMENT  = 

0. 97781 987758092D+00 

TIME 

= 

2.130, 

DISPLACEMENT  = 

0. 98242001 202845D+00 

TIME 

= 

2.140, 

DISPLACEMENT  = 

0.98701 187638821 D+00 

TIME 

= 

2.150, 

DISPLACEMENT  = 

0.991 5981 828871 4D+00 

TIME 

= 

2.160, 

DISPLACEMENT  = 

0 . 996181 57630830D+00 

TIME 

= 

2.170, 

DISPLACEMENT  = 

0.1 0007641 98791 2D+01 

TIME 

= 

2.180, 

DISPLACEMENT  = 

0. 10053494399595D+01 

TIME 

= 

2.190, 

DISPLACEMENT  = 

0. 10099393168333D+01 

TIME 

= 

2.200, 

DISPLACEMENT  = 

0. 1014536223775 0D+01 

TIME 

= 

2.210, 

DISPLACEMENT  = 

0. 10191420575858D+01 

TIME 

= 

2.220, 

DISPLACEMENT  = 

0. 102375997740670+01 

TIME 

= 

2.230, 

DISPLACEMENT  = 

0.1 028391 771 2906D+01 

TIME 

= 

2.240, 

DISPLACEMENT  = 

0. 10330396140693D+01 

TIME 

= 

2.250, 

DISPLACEMENT  = 

0.1 037705 1845379D+01 

TIME 

= 

2.260, 

DISPLACEMENT 

= 

0.1 042391 4489334D+01 

TIME 

= 

2.270, 

DISPLACEMENT 

= 

0.1 04709999323 14D+01 

TIME 

= 

2.280, 

DISPLACEMENT 

= 

0.1 05 183280291 71D+01 

TIME 

= 

2.290, 

DISPLACEMENT 

= 

0. 10565913655695D+01 

TIME 

= 

2.300, 

DISPLACEMENT 

= 

0. 10613784905069D+01 

TIME 

= 

2.310, 

DISPLACEMENT 

= 

0.1 1476631005725D+01 

TIME 

= 

2.320, 

DISPLACEMENT 

= 

0. 1 1533185695767D+01 

TIME 

= 

2.330, 

DISPLACEMENT 

= 

0.11 590500263745D+01 

TIME 

= 

2.340, 

DISPLACEMENT 

= 

0.1 1648631801910D+01 

TIME 

= 

2.350, 

DISPLACEMENT 

= 

0.11 70762 16046980+01 

TIME 

= 

2.360, 

DISPLACEMENT 

= 

0.11 7675 16974758D+01 

TIME 

= 

2.370, 

DISPLACEMENT 

= 

0.1 182836581 17380+01 

TIME 

= 

2.380, 

DISPLACEMENT 

= 

0.1 1890210681560D+01 

TIME 

= 

2.390, 

DISPLACEMENT 

= 

0.1 19531 12520959D+01 

TIME 

= 

2.400, 

DISPLACEMENT 

= 

0. 120171 15269290D+01 

TIME 

= 

2.410, 

DISPLACEMENT 

= 

0 . 1 2082269375580D+01 

TIME 

= 

2.420, 

DISPLACEMENT 

= 

0. 12148619575 064D+01 

TIME 

r 

2.430, 

DISPLACEMENT 

= 

0.1 221 62303280 190+01 

TIME 

= 

2.440, 

DISPLACEMENT 

= 

0. 12285147870224D+01 

TIME 

= 

2.450, 

DISPLACEMENT 

= 

0. 12355425415795D+01 

TIME 

= 

2.460, 

DISPLACEMENT 

= 

0.1 24271 0999 11260+01 

TIME 

= 

2.470, 

DISPLACEMENT 

= 

0.12500269925297D+01 

TIME 

= 

2.480, 

DISPLACEMENT 

= 

0. 12574953844 797D+01 

TIME 

= 

2.490, 

DISPLACEMENT 

= 

0. 12651217887352D+01 

TIME 

= 

2.500, 

DISPLACEMENT 

= 

0. 127291 11430751D+01 

TIME 

= 

2.510, 

DISPLACEMENT 

= 

0. 12808706969625D+01 

TIME 

= 

2.520, 

DISPLACEMENT 

= 

0. 12890055553986D+01 

TIME 

= 

2.530, 

DISPLACEMENT 

= 

0.1 29732 16346622D+01 

TIME 

= 

2.540, 

DISPLACEMENT 

= 

0. 13058241075455D+01 

TIME 

= 

2.550, 

DISPLACEMENT 

= 

0. 1314520664972 7D+01 

TIME 

= 

2.560, 

DISPLACEMENT 

= 

0. 13234166514682D+01 

■ 


Appendix  B 


Another  sample  program  - use  physical  units  rather  nondimensionalized  numbers. 


PROGRAM  GPLATE 

DOUBLE  PRECISION  DISPV(1024) 

DIMENSION  T I ME ( 1024) 

INTEGER  YORNO 

C 

URITEO ,*) 'THIS  PROGRAM  WILL  COMPUTE  THE  DISPLACEMENT  U SUB  I AT 
WRITEO  ,*) ' A POINT  ON  THE  SURFACE  OF  AN  INFINITE  PLATE  DUE  TO  A 
URITEO ,*) 1 POINT  STEP  IMPULSIVE  FORCE  F SUB  J OR  A POINT  LINEAR 
WRITEO  ,*) 1 RAMP  FORCE  DIPOLE  F SUB  J,K  ACTING  AT  ANOTHER  POINT1 
WRITEO,*)'  ON  THE  SURFACE  OF  THE  PLATE.1 

WRITEO,*)'  INPUTS  REQUIRED  ARE  PROMPTED  AND  TYPED  IN  ON  THE  ' 
WRITEO,*)'  KEYBOARD.' 

WRITEO,*)'  ' 

10  WRITEO,*)'  SUBSCRIPTS  OF  G = ? O.  e.  11,13,131,322,  etc.)' 
READ(1,*)  INDEX 

WRITEO,*)'  SHEAR  WAVE  SPEED  = ?' 

READ(1,*)  SHSPD 

WRITEO,*)'  LONGITUDINAL  WAVE  SPEED  = ?' 

READ(1,*)  PSPD 

100  WRITEO,*) 'ARE  SOURCE  AND  DETECTER  ON  THE  SAME  SIDE  OF  THE  PLATE 
WRITEO,*)'  PLEASE  TYPE  1 FOR  YES' 

READ  (1,*)  YORNO 
IF  (YORNO. EQ.1)  THEN 
ZD  = 0.5 
ELSE 

ZD  = -0.5 
END  IF 

WRITEO,*)'  WHAT  IS  THE  THICKNESS  OF  THE  PLATE  ?' 

READ(1,*)  HTHICK 

WRITEO,*)'  WHAT  IS  THE  DISTANCE  BETWEEN  DETECTER  AND  SOURCE' 
WRITEO,*)'  OR  EPICENTER  ?' 

READO,*)  X 

URITEO,*)'  SAMPLING  TIME  INTERVAL  = ?' 

READO,*)  DT 

WRITEO,*)'  NUMBER  OF  SAMPLING  POINTS  = ? (MAX  = 1024)' 

READO,*)  NPT 

WRITEO,*)'  SHEAR  MODULUS  OF  THE  PLATE  = ?' 

READO,*)  SHMDUL 

WRITEO,*)'  THE  OUTPUT  WILL  BE  TIME  VERSUS  DISPLACEMENT.' 
WRITEO,*)'  ' 

C 

ALPHA  = SHSPD/PSPD 
XD  = X/HTHICK 
TSCALE  = HTHICK/SHSPD 
DELTAT  = DT/TSCALE 

DSCALE  = 1 . / ( 2 . *AS I N ( 1 .0)*SHMDUL*HTHICK) 

C 

WRITE(2,*) ' ' 

URITE(2,*)'  THE  FOLLOWINGS  ARE  INPUT  PARAMETERS:' 

WRITE(2,*) ' INDEX  = '.INDEX 

WRITE(2,*) ' SHEAR  WAVE  SPEED  = ', SHSPD 

WRITE(2,*)' LONGITUDINAL  WAVE  SPEED  = ',PSPD 
WRITE (2,*) 'THICKNESS  OF  THE  PLATE  = ' , HTHICK 
WRITE(2,*) 'SAMPLING  TIME  INTERVAL  = ' ,DT 


WRITE(2,*) 1 SHEAR  MODULUS  = ' ,SHMDUL 

WRITE (2,*) 'NUMBER  OF  POINTS  TO  BE  COMPUTED  = ',NPT 
IF  (YORNO.EQ. 1 ) THEN 

WRITE(2,*)  ' THE  DETECTER  AND  THE  SOURCE  ARE  ON  THE  SAME  SIDE1 
ELSE 

WRITE(2,*)'THE  DETECTER  AND  THE  SOURCE  ARE  ON  OPPOSITE  SIDES' 
END  IF 

WRITE(2,*)  ' 1 

CALL  GREEN FCT( ALPHA, XD, ZD, INDEX,DELTAT,NPT,DISPV) 

DO  200  1=1, NPT 

TIME( I )=DELTAT*REAL( I )*TSCALE 
D I SPV( I )=DI SPV( I )*DSCALE 

WRITE(2,*) 1 TIME  = ',TIME(I),'  DISPLACEMENT  = ',DISPV(I) 
200  CONTINUE 
STOP 
END 


Appendix  C 


Listing  of  Subroutines  for  Computing  Green's  function  of  an  infinite  plate. 


c GR000020 

SUBROUTINE  GREENFCT(ALPHAf XD, ZD, INDEX, TDELTA, NPOINT, DISPL)  GR000030 

c GR000040 

q*******************-************************************************  GR000050 

c GR000060 

C TO  COMPUT  GREEN'S  FUNCTION  - DISPLACEMENT  AS  A FUNCTION  OF  TIME  GR000070 

C - FOR  GIVEN  TEST  CONFIGURATION  GR000080 

C GR000090 

C INPUTS:  GR000100 

C ALPHA  = SHEAR  WAVE  SPEED  / LONGITUDINAL  WAVE  SPEED  GR000110 

C XD  = X COORDINATE  OF  THE  DETECTOR  GR000120 

C ZD  = Z COORDINATE  OF  DETECTOR  GR000130 

C ZD  = 0.5  =>  DETECTOR  & SOURCE  ON  THE  SAME  SIDE  GR000140 

C ZD  =-0.5  =>  DETECTOR  & SOURCE  ON  OPPOSIT  SIDES  GR000150 

C INDEX  = SUBSCRIPT  OF  GREEN'S  FUNCTION,  i.e.  33,  322  GR000160 

C TDELTA  = TIME  INCREMENT  GR000170 

C NPOINT  = TOTAL  NUMBER  OF  POINTS  TO  BE  COMPUTED  GR000180 

C GR000190 

C OUTPUT:  GR000200 

C DISPL  = DISPLACEMENT  **  DOUBLE  PRECISION,  MUST  BE  GR000210 

C DIMENSIONED  (NPOINT)  GR000220 

C GR000230 

C GR000240 

C NOTE  1:  GR000250 

C ALL  INPUTS  AND  OUTPUTS  ARE  NOND I MENS I ONAL I ZED  PARAMETERS  GR000260 

C DISPL  = NORMALIZED  DISPLACEMENT  = ACTUAL  DISPLACEMENT  * PI  GR000270 

C * SHEAR  MODULUS  * PLATE  THICKNESS  / FORCE  GR000280 

C T = NORMALIZED  TIME  = ACTUAL  TIME  *SHEAR  WAVE  SPEED  / PLATE  GR000290 

C THICKNESS  GR000300 

C DISPL  MUST  BE  DIMENSIONED  IN  THE  CALLING  PROGRAM  (NPOINT)  GR000310 

C GR000320 

C GR000330 

C GR000340 

C NOTE  2:  GR000350 

C THE  ARRIVAL  TIMES  ARE  STORED  IN  RAYTIME(N, I ) GR000360 

C N = THE  Nth  RAY  GR000370 

C RAYTIME(N,3)  = THE  ARRIVAL  TIME  GR000380 

C RAYTIME(N, 1 ) = Nl,  NUMBER  OF  P TRIPS  GR000390 

C RAYTIME(N,2)  = N2,  NUMBER  OF  S TRIPS  GR000400 

C IF  N2<0,  IT  IS  A HEAD  WAVE  GR000410 

C IF  N1=- 1 , IT  IS  SURFACE  WAVE  GR000420 

C IF  N1=0  AND  N2=0,  IT  IS  A RAYLEIGH  WAVE  GR000430 

C GR000440 

C THE  TIME  OF  ARRIVAL  I M FORMAT I ON  IS  ALSO  PRINTED  ON  LU  6 GR000450 

C GR000460 

C GR000470 

C GR000480 

DOUBLE  PRECISION  DISPL(NPOINT)  GR000490 

DIMENSION  TA(101),CN(101 ,3) , RAYTIME ( 101 ,3)  GR000500 

C NRAY  = NUMBER  OF  MAXIMUM  RAYS  EXPECTED  GR000510 

DATA  NRAY/101/  GR000520 

C GR000530 

CALL  INIT(ALPHA,XD, ZD, INDEX)  GR000540 


CALL  T IMEARR I (ALPHA , XD , 2D , NRAY, RAYTIME , TA, CN ) 

WRITE(6,*)  1 I N1  N2  TIME  ARRIVAL' 

DO  500  1=1, NRAY 

WRITE  (6,499)  I , RAYTIME( I , 1 ), RAYTIME ( I , 2) , RAYTIME ( I ,3) 

499  FORMAT (2X,I5,3X,3F10.5) 

500  CONTINUE 

DC  1000  1=1 .NPOINT 

TIME  I = FLOAT ( I )*TDELTA 

CALL  GREENSUB (ALPHA , I NDEX , XD , ZD , NRAY , T I ME  I , RAYT I ME , D I SPL ( I ) ) 
WRITE(6,*) 'TIME  = ' ,T IMEI , ' DISPL  = ’,DISPL(I) 

1000  CONTINUE 
RETURN 
END 


SUBROUTINE  INIT(ALPHA,XD, ZD, INDEX) 

C 

£***★**★★*****★****★**★*★**★★*★*★*********★★*★★★***★*★******★★★**★** 


IMPLICIT  DOUBLE  COMPLEX  (Y) 

IMPLICIT  DOUBLE  PRECISION  (D) 

COMMON  /BLKO/DXD,DT,DZ,DA 

COMMON  /BLK1/KASE ,M, L , K, N1 ,N2,D,DTHETA 

COMMON  /BLK3/Y I , YONE , YTWO, YTHREE , YFOUR , YE  I GHT , YZERO 

COMMON  /BLK5/YR,Y1,Y2 

DXD  = DBLE(XD) 

DZ  = DBLE(ZD) 

DA  = DBLE(ALPHA) 

IF  (ZD.EQ.0.5)  KASE  = 1 
IF  (ZD . EQ. -0.5)  KASE  = 2 
IF  (INDEX.LT. 100)  THEN 
K = 0 

M = INDEX/10 
L = INDEX  - M*10 

ELSE  IF  ((INDEX. GT. 100). AND. (INDEX.LT. 334))  THEN 
M = INDEX/100 
L = (INDEX  - M*100)/10 
K = (INDEX  - M*100  - L*10) 

ELSE 

WRITE(6,*)  ' WRONG  INDEX  ' 

PAUSE 
END  IF 

YI  = DCMPLX(0.D0, 1 .DO) 

YONE  = DCMPLX(1 .D0,0.D0) 

YTWO  = DCMPLX(2.D0,0.D0) 

YTHREE  = DCMPLX(3.D0,0.D0) 

YFOUR  = DCMPLX(4.D0,0.D0) 

YEIGHT  = DCMPLX(8.D0, 0.D0) 

YZERO  = DCMPLX(0.D0,0.D0) 

CALL  RAYRT(YR,Y1,Y2, ALPHA) 

RETURN 

END 

C 

SUBROUT  I NE  GREENSUB( ALPHA, I NDEX , XD , ZD , NRAY , T , RAYT I ME , D I SPL ) 
C 
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£«■•*»«•»*********************•**************************************** 


C SUBROUTINE  TO  COMPUTE  THE  GREENS  FUNCTION  OF  A PLATE  AT  T 
C 

C INPUTS: 

C CONFIGURATION  PARAMETERS: 

C DETECTOR  LOCATION: 

C X = XD 

C Y = 0 

C Z = ZD,  ZD =0.5  =>  TOP;  ZD=-0.5  =>  BOTTOM 

C SOURCE  LOCATION: 

C XO  = 0 

C YO  = 0 

C ZO  = +0.5  (ON  TOP  OF  THE  PLATE) 

C 

C MATERIALS  PARAMETERS: 

C ALPHA  = RATIO  OF  SHEAR  VELOCITY  TO  LONGITUDINAL  VELOCITY 

C 

C OTHER  PARAMETERS: 

C INDEX:  SUBSCRIPTS  OF  GREENS  FUNCTION  i.e.  33,  or  113 

C NRAY : MAXIMUM  NUMBER  OF  RAYS  EXPECTED.  -100 

C RAYTIME:  TIME  OF  ARRIVAL  TABLE  OF  DIMENSION  (NRAY, 3) 

C 

C T = NORMALIZED  TIME  = ACTUAL  TIME  *SHEAR  WAVE  SPEED  / PLATE 

C THICKNESS 

C 

C DISPLACEMENT  WILL  BE  CORRECT  ONLY  IF  T < RAYTIME (NRAY, 3) 

C 

C OUTPUT: 

C DISPL  = NORMALIZED  DISPLACEMENT  = ACTUAL  DISPLACEMENT  * PI 
C * SHEAR  MODULUS  * PLATE  THICKNESS  / FORCE 

C 

C NOTE:  SUBROUTINE  INIT(ALPHA,XD, ZD, INDEX)  AND 

C SUBROUTINE  TIMEARRK ) SHOULD  BE  CALLED  FIRST  IN 

C THE  MAIN  PROGRAM. 


C 

IMPLICIT  DOUBLE  COMPLEX  (Y) 

IMPLICIT  DOUBLE  PRECISION  (D) 

DOUBLE  PRECISION  DINTG 
DIMENSION  RAYTIME(NRAY,3) 

COMMON  /BLKO/DXD,DT,DZ,DA 

COMMON  /BLK1/KASE,M,L,K,N1,N2,D,DTHETA 

COMMON  /BLK2/YD,Y, YP, YQ,YR0,YRP, YRM,YDELTA,YPSQP1 , YQSQP1 , 

1 YPSQM1 , YQSQM1 ,YSQRP1 , YSQRQ1 , YY,YA,YAA, YETA1 , YETA2 

2 , YDPHI 

COMMON  /BLK3/YI ,YONE,YTWO, YTHREE,YFOUR,YEIGHT, YZERO 
COMMON  /BLK5/YR, Y1 ,Y2 
C 

EXTERNAL  DINTG 
C 
C 

IF  ( (INDEX. EQ. 12). OR. (INDEX. EQ. 21). OR. 

+ ( INDEX. EQ. 23). OR. ( INDEX. EQ. 32). OR. 

+ ( INDEX. EQ. 121 ). OR. ( INDEX. EQ. 211). OR. 

+ ( INDEX. EQ. 231 ). OR. ( INDEX. EQ. 321). OR. 

+ ( INDEX. EQ. 123). OR. ( INDEX. EQ. 213). OR. 
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c 

c 


c 

c 

c 


c 

c 

c*** 

c 

c*** 


c*** 

c 

C8000 

8000 


( INDEX. EQ. 233). OR. ( INDEX. EO. 323). OR. 

( INDEX. EQ. 132). OR. ( INDEX. EQ. 312). OR. 

(INDEX. EQ. 112). OR. (INDEX. EQ. 222). OR. (INDEX. EQ. 332)  ) THEN 
DISPL=0.0D0 
RETURN 
ELSE 


DZERO=O.DO 

NMAX  = NRAY;  Maximum  number  of  rays.  The  same  as  DIM  RAYTIME 

NMAX  = NRAY 

IF(T.GT. (RAY TIME (NRAY , 3 ) ) ) THEN 

yp I TE  C 6 * ) * ************************************ * 
WRITE(6,*) 'ERROR  IN  SUBROUTINE  GREEN' 

WRITE(6,*)'T  > TARRIVAL(NRAY) , INCREASE  NRAY.' 
WRITE(6,*) 'PROGRAM  IS  PAUSED' 
yp I TE  C 6 * ) * *********************************** • 

PAUSE  'GREENSUB' 

END  IF 

DT=DBLE(T) 


DSURF  = 0.D0 
IF  (KASE.EQ.1)  THEN 

CALL  SFWAVE(ALPHA,XD,T, INDEX, DSURF) 

ELSE  IF  (KASE.EQ.2)THEN 
DSURF  = 0.D0 
END  IF 

DISPL=0.D0 

DO  8000  N=1 ,NMAX 

IF  (T.LE.RAYTIME(N,3))  GO  TO  8001 
N1  = RAYTIME(N , 1 ) 

N2  = RAYTIME(N, 2) 

I F ( ( N 1 . EQ.O) . AND . (N2. EQ . 0) ) GOTO  8000 
I F ( ( N 1 . LT.O) .AND . (N2.EQ.0) ) GOTO  8000 
IT  IS  RAYLEIGH  WAVE  ARRIVAL  WHICH  HAS  BEEN  TAKEN 
CARE  OF  BY  SUBROUTINE  SFWAVE 
IF(N2.LT.O)  THEN 
IT  IS  A HEAD  WAVE 
N2=-N2 

DK2=D8LE(REAL(N2) ) 

DTSARR  = DSQRT (DXD*DXD  + DK2*DK2) 

IF  (DT.GE. DTSARR)  GO  TO  8000 

DHEAD  = (DT/DA) -DK2*DSQRT ( 1 .0D0-DA*DA)/DA 

DRAYI=DINTG( -DXD, -DHEAD) 

ELSE 

IT  IS  A REGULAR  RAY 
DRAYI=DINTG(DZERO, -DXD) 

END  IF 

DISPL  = DISPL  + DRAY  I 
Y=YZERO  ;set  try  root  eq.  0 
Y=YZERO 
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8001  D I SPL=0 I SPL/(6. 2831853071 795900) 

D I SPL=DI SPL+DSURF 
RETURN 
END  IF 
END 


C 

DOUBLE  PRECISION  FUNCTION  DINTG(DINIT,DEND) 


C 


£*«***************************************************************** 


c 

IMPLICIT  DOUBLE  PRECISION  (D) 

C 

C 

DOUBLE  PRECISION  EPS, R, E,W( 50,6) , FMIN, FMAX, F 
INTEGER  NINT,NMAX,KF, I FLAG 
LOGICAL  RST 
C 
C 

EXTERNAL  F 
C 

EPS  = 1 .0-6 
RST  = .FALSE. 

NMAX  = 50 
NINT  = 1 
W(1,1)  = DINIT 
W(2, 1 ) = DEND 
C 

CALL  DGLQ1 ( F,DINIT, DEND, EPS, R,E, NINT, RST,U, NMAX, FMIN, FMAX, KF 
1 , 1 FLAG) 

C 

C 

DINTG=R 

RETURN 

END 

Q******************************************************************* 

c 

DOUBLE  PRECISION  FUNCTION  F(DS) 

C 

Q******************************************************************* 

c 

IMPLICIT  DOUBLE  COMPLEX  <Y) 

IMPLICIT  DOUBLE  PRECISION  (D) 

DOUBLE  PRECISION  ETA1 ,ETA2SQ,DD,AA 
DIMENSION  YQD(3,3),YRSTAR(3,3), YQS(3,3) 

COMMON  /BLK0/DXD,DT,DZ,DA 

COMMON  /BLK1/KASE,M,L,K,N1 ,N2,D,DTHETA 

COMMON  /BLK2/YD, Y, YP, YQ, YRO, YRP, YRM, YDELTA, YPSQP1 , YQSQP1 , 

1 YPSQM1 , YQSQM1 , YSQRP1 , YSQRQ1 , YY, YA , YAA , YETA1 , YETA2 

2 , YDPHI 

COMMON  /BLK3/Y I , YONE , YTWO, YTHREE , YFOUR , YE  I GHT , YZERO 
C 

DK1=DBLE(REAL(N1 ) ) 

DK2=DBLE(REAL(N2) ) 

D=DS 


GR002200 

GR002210 

GR002220 

GR002230 

GR002240 

GR002250 

GR002260 

GR002270 

GR002280 

GR002290 

GR002300 

GR002310 

GR002320 

GR002330 

GR002340 

GR002350 

GR002360 

GR002370 

GR002380 

GR002390 

GR002400 

GR002410 

GR002420 

GR002430 

GR002440 

GR002450 

GR002460 

GR002470 

GR002480 

GR002490 

GR002500 

GR002510 

GR002520 

GR002530 

GR002540 

GR002550 

GR002560 

GR002570 

GR002580 

GR002590 

GR002600 

GR002610 

GR002620 

GR002630 

GR002640 

GR002650 

GR002660 

GR002670 

GR002680 

GR002690 

GR002700 

GR002710 

GR002720 

GR002730 

GR002740 


c 

CALL  PHIEQ0(DK1 , DK2, DA, DT,DXD,DS,Y, DERR, NROOT) 

C WRITE(6,*)DK1  ,DK2,DS,Y, DERR, NROOT 

IF  (NROOT. EQ.O)  THEN 
F=0.D0 
C 

ELSE 

CALL  COMMON (DS) 

CALL  RSTAR(YRSTAR, IDUMMY) 

CALL  QD(YQD, IDUMMY) 

CALL  QSU(YQS, IDUMMY) 

C 

YINTG  = YZERO 
DO  200  1=1,2 
DO  100  J=1 , 2 

YINTG=YINTG+YQD(M, I )*YRSTAR( I , J )*YQS(L, J ) 
100  CONTINUE 

200  CONTINUE 

C 

IF  (<N1. EQ.O). AND. (N2.GT.0)) 

1 YINTG=YINTG+YQD(M,3)*YQS(L ,3) 

C 

YINTG=YINTG*YDPHI/CDSQRT (DCMPLX(DXD*DXD-DS*DS) ) 
F=DBLE(YINTG) 

C 

END  IF 
RETURN 
END 


C 

SUBROUTINE  PHIEQ0(DK1 ,DK2,DA,DT,DXD, D,YROOT, DERR, NROOT) 
C 


C 

C GIVEN  DK1 ,DK2,DA,DT,D  TO  SOLVE  PHI(YROOT)=0 
IMPLICIT  DOUBLE  PRECISION  (D) 

IMPLICIT  DOUBLE  COMPLEX  (Y) 

DOUBLE  PRECISION  T,R,A,C1,C2 
COMMON  /PARM1/T,R,A,C1 ,C2 
C 

YOLD=YROOT 

IER=0 

DYREAL=O.DO 

YPHI=DCMPLX(0.D0,0.D0) 

T = DT 
R = D 
A = DA 
Cl  = OKI 
C2  = DK2 
EPS=1 . E - 1 2 
IEND=100 
C 

NROOT=1 

IF  (DK2.EQ.0.D0)  THEN 
C 
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DSTSQ  = (DT/DA)*(DT/DA) -DK1*DK1 
DSTAR  = DSQRT (DSTSQ) 

IF  ( D.LT.O. DO )DSTAR=- DSTAR 

IF  (D.EQ.O.DO)  THEN 

DYIMAG=-DK1/(DSTAR*DA) 

YROOT  = DCMPLX(O.DO,DYIMAG) 

ELSE  IF  (D.EQ. DSTAR)  THEN 
DYREAL=DT/(DA*DA*DSTAR) 

YROOT  = DCMPLX(DYREALf O.ODO) 

ELSE 

DYREAL  = D*DT/(DA*DA*DSTSQ) 

DSMD=DSTSQ-D*D 

IF(DSMD.LT.O.DO)  WRITE(6,*)  DT,D, DSHD, DSTSQ 
DYIMAG  = * (DK1/(DSTSQ*DA)  )*DSQRT (DSMD) 

YROOT  = DCMPLX(DYREAL, DYIMAG) 

END  IF 

YOLD  = YROOT 

CALL  DCRTNI (YROOT, YPHI , YDPHI , YOLD, EPS, I END, IER) 

LSE  IF  (DK1.EQ.0.D0)  THEN 

DTHARR  = DA*DABS(DXD)+DK2*DSQRT(1.0D0-DA*DA) 
DTSARR  = DSQRT (DXD*DXD+DK2*DK2) 

DSTSQ  = (DT)*(DT) -DK2*DK2 
DSTAR  = DSQRT (DSTSQ) 

IF  (DT.GT. DTSARR)  THEN 

IF  (D.LT.O. DO)  DSTAR=-DSTAR 

IF  (D.EQ.O.DO)  THEN 

DYIMAG=-DK2/DABS(DSTAR) 

YROOT  = DCMPLX(0. DO, DYIMAG) 

ELSE  IF  (D.EQ. DSTAR)  THEN 
DYREAL=OT/ (DSTAR) 

YROOT  = DCMPLX(DYREAL, O.ODO) 

ELSE 

DYREAL  = D*DT/(DSTSQ) 

DYIMAG  = - (DK2/DSTSQ)*DSQRT (DSTSQ-D*D ) 

YROOT  = DCMPLX(DYREAL, DYIMAG) 

END  IF 

ELSE  IF  ((DT.LT. DTSARR). AND. (DT.GT. DTHARR))  THEN 
DHEAD  = (DT/DA)-DK2*DSQRT(1.0DO-DA*DA)/DA 
IF  (D. LT. 0. DO )DHEAD=- DHEAD 
IF  (D.LT.O.DO)DSTAR=-DSTAR 

IF  (D.EQ.O.DO)  THEN 
DYIMAG=-DK2/DSTAR 
YROOT  = DCMPLX(0. DO, DYIMAG) 

ELSE  IF  (D.EQ. DSTAR)  THEN 
DYREAL=DT/(DSTAR) 

YROOT  = DCMPLX(DYREAL, O.ODO) 

ELSE  IF  (DABS (D ) .GT. DABS (DSTAR ) ) THEN 
IF  (D.LT.O. DO)  THEN 

DYREAL  = (D*DT-DK2*DSQRT (D*D - DSTSQ ))/( DSTSQ) 
ELSE  IF  (D.GT.O.DO)  THEN 
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DYREAL  = (D*DT+DK2*DSQRT(D*D-DSTSQ))/(DSTSQ) 
ENDIF 

DYIMAG  = O.DO 

YROOT  = DCMPLX<DYREAL, DYIMAG) 

END  IF 
END  IF 

YOLD  = YROOT 

CALL  DCRTNI (YROOT, YPHI , YDPHI , YOLD, EPS, IEND, IER) 

C 

ELSE  IF  ((DK1. NE. O.DO). AND. (DK2.NE. O.DO))  THEN 

IF((D.EQ.O.DO).OR.(YOLD.EQ.DCMPLX(O.DO,O.DO)))  THEN 
CALL  DEQO(DT, DA, DK1,DK2, YROOT, NOROOT) 

ENDIF 

YOLD  = YROOT 

CALL  DCRTNI (YROOT, YPH I , YDPHI , YOLD, EPS, IEND, IER) 

END  IF 

I F( IER.NE.O)  NROOT=0 
C 

C TEST  ROOT  TO  FIND  ERROR 
C 

DERR=CDABS(YPHI ) 

C 


RETURN 

END 


C 

SUBROUTINE  DEQO (DT , DA , DK 1 , DK2 , Y , NOROOT ) 

C 

£******************************************************************* 


C 

C COMPUTE  THE  ROOT  BY  SOLVING  QUARTIC  EQ. 
C 


IMPLICIT  DOUBLE  PRECISION  (D) 

IMPLICIT  DOUBLE  COMPLEX  (Y) 

DOUBLE  PRECISION  B 

DIMENSION  DB(5) ,B(3),DROOT1 (2) ,DROOT2(2) 

C 

YF( Y)=-DCMPLX(DT , O.DO)*Y-DCMPLX(DK1 ,O.DO)* 

1 CDSQRT (Y*Y*DCMPLX(DA*DA,O.DO) -DCMPLX( 1 .DO, O.DO) ) 

2 -DCMPLX(DK2,0.D0)*CDSQRT (Y*Y-DCMPLX( 1 .DO, O.DO) ) 
DT2  = DT*DT 

DK12  =DK1*DK1 
DK22  = DK2*DK2 
DAA=DA*DA 
C 

DB( 1 )=( (DT2-DK12*DAA-DK22)**2  - 4.D0*DK12*DK22*DAA) 
DB(3)=2.D0*((DT2-DK12*DAA-DK22)*(DK12+DK22) 

1 +2.D0*DK12*DK22*(DAA+1.D0)) 

DB(5 )=(DK12+DK22)**2-4.D0*DK12*DK22 
C 

B ( 1 ) = DB( 1 ) 

B(2)  = DB(3) 

B C 3 ) = DB(5) 

T0L=1.E-20 

C CALL  QUADRTIC  EQ.  TO  FIND  Y**2 
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CALL  QUADRT ( B , DROOT 1 , DROOT2 , TOL , NOROOT ) 
C 

Yl=CDSORT (DCMPLX(DROOTl ( 1 ), DROOT 1 (2) ) ) 
Y2=CDSQRT (DCMPLX(DROOT2( 1 ) ,DROOT2(2) ) ) 

IF  (DIMAG(YI).GT.O.DO)  Y1=-Y1 
IF  (DIMAG(Y2) .GT.O.DO)  Y2=-Y2 
YF1=YF(Y1 ) 

YF2=YF(Y2) 

DERR1=CDABS(YF(Y1 )) 

DERR2=CDABS(YF(Y2)) 

IF  (DERR1.LT. DERR2)  Y=Y1 
IF  (DERR1 .GE.DERR2)  Y=Y2 
C 

IF  (DIMAG(Y) .GE.O.DO)  Y=DCONJG(Y) 

C 


RETURN 

END 


C 

DOUBLE  COMPLEX  FUNCTION  PHIO(Z) 


C 


c 

C FUNCTION  PHIO 
C 

DOUBLE  PRECISION  T,RfA,C1,C2 
DOUBLE  COMPLEX  Z, YAA, YONE 
COMMON  /PARM1/T,R,A,C1 ,C2 
C 

YONE=DCMPLX( 1 .DO,O.DO) 

YAA=DCMPLX(A*A, 0.D0) 

C 

PHIO=DCMPLX(R, 0.D0) -DCMPLX(T,0.D0)*Z 

1 *DCMPLX(C1 , O.DO)*CDSQRT (YAA*Z*Z- YONE) 

2 -DCMPLX(C2,O.DO)*CDSQRT(Z*Z-YONE) 

C 


RETURN 

END 

£***************************★*************-************************** 


c 

DOUBLE  COMPLEX  FUNCTION  DPHIO(Z) 
C 


c 

C FUNCTION  DPHIO 
C 

c 

DOUBLE  PRECISION  T,R,AfCl,C2 
DOUBLE  COMPLEX  Zf YAA, YONE 
COMMON  /PARM1/T,R,A,C1,C2 
C 

YONE=DCMPLX(1.DO,O.DO) 
YAA=DCMPLX(A*A, 0.D0) 

C 
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DPHI0=-DCMPLX(T,0.D0) 


1 - DCMPLX ( Cl , 0 . DO )*YAA*Z/CDSQRT ( YAA*Z*Z - YONE ) 

2 -DCMPLX(C2,O.DO)*Z/CDSQRT(Z*Z-YONE) 

C 

RETURN 

END 

Q******************************************************************* 

c 

SUBROUTINE  COMMON(DS) 

C 

C 

C TO  COMPUTE  PARAMETERS  DEPENDING  ON  Y AND  TO  PUT  IN 
C A COMMON  BLK2 
C 

IMPLICIT  DOUBLE  COMPLEX  (Y) 

IMPLICIT  DOUBLE  PRECISION  (D) 

DOUBLE  PRECISION  ETA1 ,ETA2SQ,DD,AA 

COMMON  /BLKO/DXD,DT,DZ,DA 

COMMON  /BLK1/KASE,M,L,K,N1 ,N2,D,DTHETA 

COMMON  /BLK2/YD, Y, YP, YQ,YR0, YRP, YRM,YDELTA, YPSQP1 , YQSQP1 , 

1 YPSQM1 , YQSQM1 , YSQRP1 .YSQRQl , YY,YA, YAA,YETA1 , YETA2 

2 , YDPHI 

COMMON  /BLK3/YI ,YONE, YTUO, YTHREE,YFOUR, YEIGHT,YZERO 
C 

AA=DA*DA 

DD=DS 

D=DS 

DN1=DBLE(REAL(N1 ) ) 

DN2=DBLE(REAL(N2) ) 

YZERO=DCMPLX(O.DO,O.DO) 

YI=DCMPLX(0.D0, 1 .ODO) 

YAA=DCMPLX(AA,O.DO) 

YA=DCMPLX(DA,0.D0) 

YONE=DCMPLX(1 .ODO.O.ODO) 

YTUO=OCMPLX( 2.00,0. DO) 

YFOUR=DCMPLX(4.0D0,0.D0) 

C 

ETA1=DS/DXD 
ETA2SQ=1.D0-ETA1*ETA1 
YETA1  = DCMPLX(ETA1 ,0.D0) 

YETA2S=  DCMPLX (ETA2SQ,0. DO) 

YETA2=CDSQRT(YETA2S) 

C 

YD  = DCMPLX(DD,0.D0) 

YY=Y*Y 

YP=CDSQRT ( YAA*YY • YONE ) 

YQ=CDSQRT(YY- YONE) 

YPSQP1=YAA*YY 
YQSQP1=YY 
YQSQM1=YQ*YQ- YONE 
YPSQM1 =YP*YP - YONE 
C 

YSQRP1=Y*YA 

YSQRQ1=Y 

YDELTA=YQSQM1*YQSQM1+YFCUR*YP*YQ 
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YM=(YQSQM1*YQSQM1 - YFOUR*YP*YQ) 
YRO=YM/YDELTA 

YRP=- (YF0UR*YP*YQSQM1 )/( YA*YDELTA) 
YRM=( YR0*YR0- YONE)/YRP 


YDPHI=Y0NE/(YD+(DCMPLX(DN1 , 0.D0)/YP)+(DCMPLX(DN2,0.D0)/YQ) ) 

RETURN 

END 


C 


SUBRCXIT I NE  RSTARCYRSTAR, IDUMMY) 


C 


c 

C TO  COMPUTE  RSTAR  (JS,J2)  EQ.6.39 

C 


IMPLICIT  DOUBLE  COMPLEX  (Y) 

IMPLICIT  DOUBLE  PRECISION  (D) 

DIMENSION  YRSTAR(3,3),KSUM(20) 

COMMON  /BLK0/DXD,DT,DZ,DA 

COMMON  /BLK1/KASE,M,L,K, N1 ,N2,D,DTHETA 

COMMON  /BLK2/YD, Y, YP, YQ, YRO, YRP,YRM,YDELTA, YPSQP1 , YQSQP1 , 

1 YPSQM1 , YQSQM1 , YSQRP1 ,YSQRQ1 ,YY,YA, YAA, YETA1 , YETA2 

2 , YDPHI 

COMMON  /BLK3/YI , YONE, YTWO, YTHREE, YFOUR, YEIGHT,YZERO 


EXTERNAL  KDELTA 


C 

C 

IDUMMY=0 

YRSTAR(3,3)  = YONE 
YRSTAR(1 ,3)  = YZERO 
YRSTAR(2,3)  = YZERO 
YRSTAR(3, 1 ) = YZERO 
YRSTAR(3,2)  = YZERO 
C 

DO  2000  JS=1 , 2 
DO  1000  J2=1 ,2 

IF  (J2.E0.1)  N1=N1 - 1 
IF  ( J2.EQ.2)  N2=N2-1 
ISIGN  = <-1)**(J2*(N1+N2)+N2) 

CALL  COEFF(N1,N2, JS, J2,KDf KU,KL,KSUM) 

C WRITE(7,*)  N1 , N2, JS, J2,KD, KU,XL,KSUM 

YSUM=YZERO 
DO  100  NM=<L , KU 

YSUM=YSUM+DCMPLX<DBLE(REAL(KSUM(NM+1 ) ) ) , 0.D0)* 
1 ( (YR0*YR0- YONE)/(YRO*YRO) )**NM 

100  CONTINUE 

KD2=KDELTA( 1,JS)-KDELTA(1,J2) 

C 

YSUM  = ( (YR0/YRP)**(KD2)  )*YSUM 
RSI GN=REAL( ISIGN) 

DRSIGN=RSIGN 

DKD=DBLE(REALOCD)) 

YRSTAR( JS, J2)=DCMPLX(DRSI GN, 0.D0)*(YR0**(N1+N2) ) 
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1 *(DCMPLX(DKD,O.DO)+YSUM) 

IF  (J2.EQ.1)  N1=N1+1 
IF  (J2.EQ.2)  N2=N2+1 
C 

1000  CONTINUE 
2000  CONTINUE 
RETURN 
END 

c******************************************************************* 

c 

SUBROUTINE  QD(YQD, IDUMMY) 

C 

C 

C TO  COMPUTE  QD( I ( J ) BASED  ON  EQ.  6.21 
C 

IMPLICIT  DOUBLE  COMPLEX  (Y) 

IMPLICIT  DOUBLE  PRECISION  (D) 

COMMON  /BLKO/DXD,DT,DZ,DA 

COMMON  /BLK1/KASE,M,L,K,N1 ,N2,D,DTHETA 

COMMON  /BLK2/YD, Y,YP,YQ,YRO,YRP, YRM, YDELTA, YPSQP1 ,YQSQP1 , 

1 YPSQM1 , YQSQM1 , YSQRP1 , YSQRQ1 , YY , YA , YAA , YETA1 , YETA2 

2 , YDPHI 

COMMON  /BLK3/Y I , YONE , YTUO, YTHREE , YFOUR , YE  I GHT , YZERO 
C 

DIMENSION  YQD(3,3) 

C 

1DUMMY=0 
DO  100  1=1,3 
DO  90  J=1 ,3 

YGD(I,J)  = YZERO 
90  CONTINUE 

100  CONTINUE 
C 

IF(M.EQ.I)  THEN 

YQD< 1 , 1 )=YETA1*YFOUR*YP*YQ*YQSQP1 
1 / (YSQRP1*YDELTA) 

YQD( 1 ,2)=- YTW0*YETA1*YQ*YQSQM1*YSQRQ1 
1 /(YDELTA) 

YQD( 1,3)  =-  YTWO*YETA2 
C 

ELSE  IF  (M.EQ.2)THEN 
YQD(2, 1 )=YETA2*YF0UR*YP*YQ*YQSQP1 
1 /( YSQRP1*YDELTA) 

YQD(2,2)=-YTW0*YETA2*YQ*YQSQM1*YSQRQ1 
1 /(YDELTA) 

YQD(2,3)=YTW0*YETA1 

C 

ELSE  IF  (M.EQ.3)  THEN 
YQD(3 , 1 )=- YTU0*YP*YQSQM1*YQSQP1 
1 /( YSQRP1*YDELTA) 

YQD(3,2)=-YF0UR*YP*YQ*YSQRQ1 
1 /(YDELTA) 

YQD(3,3)=YZERO 
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END  IF 
RETURN 
END 

c******************************************************************* 

C 

INTEGER  FUNCTION  KDELTAd , J) 

C 

C******************************************************************* 

c 

IF  (I.EQ.J)  THEN 
KDELTA=1 
ELSE 

KDELTA=0 
END  IF 
RETURN 
END 

Q******************************************************************* 

c 

INTEGER  FUNCTION  IBINO(M,N) 

C 


IF  (M.EQ.N)GO  TO  11 

IF  ( (M. LT .0) .OR . (N. LT.O) ) GO  TO  10 

K=M-N 

IF  (1C)  10,11,12 

10  IBINO  = 0 
RETURN 

11  IBINO  = 1 
RETURN 

12  IF  (N.EQ.O)GO  TO  11 
IN  = 1 

IP  = M 
IQ  = 1 
DO  20  1=1, < 

IN  = IN* IP/ IQ 
IP  = IP-1 
IQ  = IQ  +1 
20  CONTINUE 

IBINO  = IN 
RETURN 
END 


INTEGER  FUNCTION  KD1 (N1 ,N2, JS, J2) 


Q******************************************************************* 

C 

EXTERNAL  KDELTA 

KD1=KDELTA<0,N1)*KDELTA(2,J2)+KDELTA<0,N2)*KDELTA(1,J2) 

KD1=ICD1*KDELTA<JS,J2) 

RETURN 

END 

c******************************************************************* 

c 
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INTEGER  FUNCTION  KUPPER(N1 ,N2, JS, J2) 
C 


EXTERNAL  KDELTA 

K=N2+KDELTA(2, J2)*KDELTA( 1 , JS) -KDELTA( 1 , J2)*KDELTA(2, JS) 
IF  (Nl.LE.K)  THEN 
KUPPER=N1 
ELSE 

KUPPER=K 
END  IF 
RETURN 
END 


INTEGER  FUNCTION  KLOWER(N1 , N2, JS , J2) 


EXTERNAL  KDELTA 

KLOUER=KDELTA(2, J2)+KDELTA( 1 , J2)*KDELTA( 1 , JS) 

RETURN 

END 

£★ ★★**★★***★★*★★*★*★*★★★*★*★★★ ************************************** 


C 

SUBROUTINE  COEFF(Nl , N2, JS, J2,KD,KU,KL,KSUM) 

C 

C 

C TO  COMPUTE  THE  COEFFICIENTS  FOR  THE  RSTAR  TURM 

C 

C INPUTS: 

C N1  = NUMBER  OF  P TRIPS 

C N2  = NUMBER  OF  S TRIPS 

C JS  = FINAL  TRIP  TYPE,  1 OR  2 

C J2  = INITIAL  TRIP  TYPE,  1 OR  2 

C 

C OUTPUT: 

C KD  = FIRST  COEFF 

C KU  = UPPER  LIMIT  IN  SUMMATION 

C KL  = LOWER  LIMIT  IN  SUMMATION 

C KSUM  = COEFF.  IN  SUMMATION,  MUST  BE  DIMENSIONED  IN  CALLING 

C PROGRAM  TO  BE  MAX(N1,N2) 

C 

C 

c 

C PROGRAM  TESTCOEF 

C DIMENSION  KSUM(20) 

cc*** 

C 1 WRITE(5, 10) 

C 10  FORMAT  (2X, 1 TYPE  K1 ,K2, JS, J2 • ) 

C READ(5,*)I,J,JS,J2 

C WRITE(6,30)I,J,JS,J2 

C 30  FORMAT (2X, 1 I =• , 15, 1 J= ■ , 15 , • JS=',I2,'  J2=',I2) 

CC*** 
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END  IF 
RETURN 
END 

q* a***************************************************************** 

c 

INTEGER  FUNCTION  KDELTA(I.J) 

C 


IF  (I.EQ.J)  THEN 
KDELTA=1 
ELSE 

KDELTA=0 
END  IF 
RETURN 
END 


INTEGER  FUNCTION  IBINO(H,N) 


IF  (M.EQ.N)GO  TO  11 

IF  ((M.LT.O).OR.(N.LT.O))  GO  TO  10 

K=M-N 

IF  (K)  10,11,12 

10  IBINO  = 0 
RETURN 

11  IBINO  = 1 
RETURN 

12  IF  (N.EQ.O)GO  TO  11 
IN  = 1 

IP  = M 
IQ  = 1 
DO  20  1=1, K 
IN  = IN* IP/ IQ 
IP  = IP-1 
IQ  = IQ  +1 
20  CONTINUE 

IBINO  = IN 
RETURN 
END 

£★★*★★★**★★★★★**★**★*★*★**★*★*****★**★**★★★★★****★★**★*★★★★*★★★★*★*★ 

C 

INTEGER  FUNCTION  KD1 (N1 , N2, JS, J2) 

C 


EXTERNAL  KDELTA 

KD1=KDELTA(0,N1)*KDELTA(2,J2)+KDELTA(0,N2)*KDELTA(1,J2) 
KD1=KD1*KDELTA( JS, J2) 

RETURN 

END 

£*******************************************-********* *************** 
C 
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INTEGER  FUNCTION  KUPPER(N1 , N2, JS, J2) 

C 

C******************************************************************* 

C 

EXTERNAL  KDELTA 

K=N2+KDELTA(2, J2)*KDELTA( 1 , JS)-KDELTA( 1 , J2)*KDELTA(2, JS) 

IF  (Nl.LE.K)  THEN 
KUPPER=N1 
ELSE 

KUPPER=K 
END  IF 
RETURN 
END 

c 

INTEGER  FUNCTION  KLOWER(N1 ,N2, JS, J2) 

C 

Q******************************************************************* 

c 

EXTERNAL  KDELTA 

KLOUER=KDELTA<2, J2)+KDELTA(1,J2)*KDELTA(1, JS) 

RETURN 

END 

c******************************************************************* 

c 


c 

SUBROUTINE  COEFF(N1 ,N2, JS, J2,KD,KU, 

KL,KSUM) 

c 

c 

TO  COMPUTE  THE  COEFFICIENTS  FOR  THE 

RSTAR  TURM 

c 

c 

INPUTS: 

c 

N1  = NUMBER  OF  P TRIPS 

c 

N2  = NUMBER  OF  S TRIPS 

c 

JS  = FINAL  TRIP  TYPE,  1 OR  2 

c 

J2  = INITIAL  TRIP  TYPE,  1 OR  2 

c 

c 

OUTPUT: 

c 

KD  = FIRST  COEFF 

c 

KU  = UPPER  LIMIT  IN  SUMMATION 

c 

KL  = LOWER  LIMIT  IN  SUMMATION 

c 

KSUM  = COEFF.  IN  SUMMATION,  MUST 

BE  DIMENSIONED  IN  CALLING 

c 

PROGRAM  TO  BE  MAX(N1,N2) 

c 

r... 

c 

C PROGRAM  TESTCOEF 

C DIMENSION  KSUM(20) 

CC*** 

C 1 WRITE(5, 10) 

C 10  FORMAT  (2X, 1 TYPE  K1 ,K2, JS, J2' ) 

C READ (5,*) I , J, JS, J2 

C URITE(6,30)I , J, JS, J2 

C 30  FORMAT (2X, 1 I =• , 15, ■ J=* , 15, ■ JS=',I2,'  J2=',I2) 
CC*** 


GR007150 

GR007160 

GR007170 

GR007180 

GR007190 

GR007200 

GR007210 

GR007220 

GR007230 

GR007240 

GR007250 

GR007260 

GR007270 

GR007280 

GR007290 

GR007300 

GR007310 

GR007320 

GR007330 

GR007340 

GR007350 

GR007360 

GR007370 

GR007380 

GR007390 

GR007400 

GR007410 

GR007420 

GR007430 

GR007440 

GR007450 

GR007460 

GR007470 

GR007480 

GR007490 

GR007500 

GR007510 

GR007520 

GR007530 

GR007540 

GR007550 

GR007560 

GR007570 

GR007580 

GR007590 

GR007600 

GR007610 

GR007620 

GR007630 

GR007640 

GR007650 

GR007660 

GR007670 

GR007680 

GR007690 


IF  (KASE.EQ.1)  THEN 
C SOURCE  AT  TOP  SURFACE 

C 

DO  200  1=1,3 

YQS( I , 1 )=YRRP( 1 , 1 )*YUP( I , 1 )*YP1+YRRP( 1 ,2)*YUP( I , 2)*YP2 
1 - YUM( I , 1 )*YP1 

YQS( I ,2)=YRRP(2, 1 )*YUP( I , 1 )*YP1 
1 +YRRP(2,2)*YUP(I,2)*YP2  - YUM(I,2)*YP2 

YQS(I,3)=(YRRP(3,3)*YUP(I,3)  - YUM( I ,3))*YP2 
200  CONTINUE 

C 
C 

ELSE  IF  (KASE.EQ.2)  THEN 
C SOURCE  AT  BOTTOM  SURFACE 

C 

DO  300  1=1,3 

YQS( I , 1 )=YRRM( 1 , 1 )*YUM( I , 1 )*YP1+YRRM< 1 ,2)*YUM( I , 2)*YP2 
1 - YUP( I , 1 )*YP1 

YQS( I , 2)=YRRM(2, 1 )*YUM( I , 1 )*YP1 
1 +YRRM(2,2)*YUM( I ,2)*YP2  - YUP<I,2)*YP2 

YQS( I ,3)=(YRRM(3,3)*YUM( 1,3)  - YUP(I,3)  )*YP2 
300  CONTINUE 

C 

END  IF 

RETURN 

END 

c ★★*★★★★**★★★****★*★**★★***★★****★******★★*******★*★**★★********★★ 

C 

SUBROUTINE  EPIDIS(A,T,U) 

C 

Q *★★★★*★★*★*★★★★*★★★★**★*★★★★*★*★★*★★★★★★*★*****★**★★★****★*★★★**★ 

C 

C 

C THIS  PROGRAM  CALCULATES  THE  VERTICAL  DISPLACEMENT 

C AT  THE  EPICENTER  OF  A LARGE  PLATE  DUE  TO  AN  UNIT 

C STEP  FUNCTION  VERTICAL  FORCE. 

C DOUBLE  PRECISION  IS  USED. 

C 

C INPUT  PARAMETERS  : 

C 

C A = SHEAR  WAVE  SPEED  / LONGITUDINAL  WAVE  SPEED 
C e.g.  ALUMINUM  : A=0.495 

C T = NORMALIZED  TIME  = ACTUAL  TIME  *SHEAR  WAVE  SPEED  / PLATE 
C THICKNESS 

C e.g.  SHEAR  WAVE  SPEED  ~ 3 mm  /micro  second 
C 

C OUTPUT  PARAMETERS  : 

C 

C U = NORMALIZED  DISPLACEMENT  = ACTUAL  DISPLACEMENT  * 3.14159 
C * SHEAR  MODULUS  * PLATE  THICKNESS  / FORCE 
C 

C SUBROUTINES  REQUIRED: 

C 

C QPARA(T,A>,  AND  FUNCTION  BINO 
C 


GR008800 
GR008S10 
GR008820 
GR008830 
GR008840 
GR008850 
GR008860 
GR008870 
GR008880 
GR008890 
GR008900 
GR008910 
GR008920 
GR008930 
GR008940 
GR008950 
GR008960 
GR008970 
GR008980 
GR008990 
GR009000 
GR009010 
GR009020 
GR009030 
GR009040 
GR009050 
GR009060 
GR009070 
GR009080 
GR009090 
GR009100 
GR0091 10 
GR009120 
GR009130 
GR009140 
GR009150 
GR009160 
GR009170 
GR009180 
GR009190 
GR009200 
GR009210 
GR009220 
GR009230 
GR009240 
GR009250 
GR009260 
GR009270 
GR009280 
GR009290 
GR009300 
GR009310 
GR009320 
GR009330 
GR009340 


c- 

GR009350 

c 

GR009360 

c 

THE  PROGRAM  MAY  BE  SET  UP  AS  THE  FOLLOWINGS: 

GR009370 

c 

GR009380 

c- 

PROGRAM  EPICT 

GR009390 

c- 

DOUBLE  PRECISION  A,DT,DELTAT,DU(1024) 

GR009400 

c- 

A = 0.584D0 

GR009410 

c- 

NPT=1024 

GR009420 

c- 

DELTAT=0.00277D0 

GR009430 

c- 

DO  100  1=1 f NPT 

GR009440 

c- 

DT=DELTAT*DBLE(REAL( I ) ) 

GR009450 

c- 

CALL  EPIDIS(A,DT,DU( I ) ) 

GR009460 

c- 

100  CONTINUE 

GR009470 

c- 

WRITE(7)  DU 

GR009480 

c- 

STOP 

GR009490 

c- 

END 

GR009500 

c 

GR009510 

GR009520 

c 

GR009530 

IMPLICIT  DOUBLE  PRECISION  (Q) 

GR009540 

DOUBLE  PRECISION  U,T,A,TA,DXM 

GR009550 

INTEGER  BINO 

GR009560 

COMMON  /BLOCK/  GN1 ,QN2,QO,QP,QM,QPP,QPS,QSP,QSS,QPHI 

GR009570 

c 

GR009580 

c 

NRAY  = MAX.  POSSIBLE  NUMBER  OF  P RAYS 

GR009590 

c 

NRAY2  = MAX.  POSSIBLE  NUMBER  OF  S RAYS 

GR009600 

NRAY  = 20 

GR009610 

NRAY2  = 12 

GR009620 

U = 0.0 

GR009630 

DO  300  NA  = 1 , NRAY 

GR009640 

DO  200  NB  = 1.NRAY2 

GR009650 

N1  = NA- 1 

GR009660 

N2  = NB- 1 

GR009670 

c 

N1+N2  MUST  BE  AN  ODD  NUMBER 

GR009680 

IF  ( ( FLOAT (N1+N2)/2. ) .EQ. FLOAT ( (N1+N2)/2) )GO  TO  200 

GR009690 

TA  = N2  + N1*A 

GR009700 

IF  (T.LE.TA)  GO  TO  200 

GR009710 

QN1  = N1 

GR009720 

QN2  = N2 

GR009730 

c 

GR009740 

CALL  QPARA(T,A) 

GR009750 

c 

GR009760 

NTERM  = N2 

GR009770 

IF  (N1.LT.N2)  NTERM=N1 

GR009780 

NTERM=NTERM+1 

GR009790 

QSUM  = 0.0 

GR009800 

DO  100  MT  = 1, NTERM 

GR009810 

MM  = MT- 1 

GR009820 

QFPP  = BINO(N1+1,MM+2)*BINO(N2+1,MM+1) 

GR009830 

QFPS  = BIN0(N2+1,MM+1)*BIN0(N1+1,MM+1) 

GR009840 

QFSS  = BINO(N2+1 ,MM+2)*BIN0(N1+1 ,MM+1 ) 

GR009850 

MO  = N1  + N2  - 1 - 2*MM 

GR009860 

c 

GR009870 

QSUM  = ( ( - 1 . )**N2)* 

GR009880 

1 ( (QFPP*QPP-QFSS*QSS)*(DXM(QO, MO)*DXM( (QM*QP*( - 1 . ) ) , MM) ) 

GR009890 

2 +QFPS*(QPS*DXM(QO, (MO+1 ) )*DXM( ( ( - 1 . )*QP) ,MM)*DXM(QM,MM) 

3 -OSP*DXM(QO, (MO+1 ) )*OXM( ( ( - 1 . )*QM) ,MM)*DXM(QP, MM) ) ) 

4 +QSUM 
C 

100  CONTINUE 

U = U + QSUM/QPHI 
200  CONTINUE 
300  CONTINUE 
RETURN 
END 
C 

Q ************************** *************************************** 

c 

SUBROUTINE  QPARA(T.A) 

C 

Q ★★★★★★★*★★★**★★*★★★★★★★**★*★★*★****★★★★*★*★★*★★★★★★★****★****★★★★ 

C 

IMPLICIT  DOUBLE  PRECISION  (Q) 

DOUBLE  PRECISION  T^B.AC.BAC.XX.X.SX.SAX 
COMMON  /BLOCK/  QN1 , QN2,QO, QP,QM,QPP,QPS,QSP,QSSf QPHI 
C 

C CALCULATE  X 
C 

Q1=QN1*QN1 
Q2=QN2*QN2 
IF(QNI.EQ.O.DO)  THEN 
XX=(T*T-Q2)/Q2 
ELSE  IF  (QN2.EQ.0.D0)  THEN 
XX=(T*T-Q1*A*A)/Q1 
ELSE 

B = ((Q1+Q2)*T*T  -(Q1*A*A-Q2)*(Q1-Q2))/((Q1-Q2)**2) 

AC=  ( (Q1*A*A-Q2)**2+T*T*(T*T - 2 . *Q1 *A*A- 2 . *Q2) ) 

1 /( (Q1 -Q2)**2) 

BAC  =B*B  -AC 

IF  (BAC. LE. 0.0)  3AC  = 0.0 
XX  = B-DSQRT (BAC) 

END  IF 

X = DSQRT(XX) 

SX  = DSQRT ( 1 .+XX) 

SAX  = DSQRT (A*A+XX) 

C 

C CALCULATE  R'S,  Q'S,  AND  QPHI 
C 

QA=  (1.+2.*XX)**2 
QB=  4.*SAX*SX*XX 
QO=  (QA+QB)/(QA-QB) 

QP=  -4.*X*A*SX*(1.+2.*XX)/(QA-QB) 

QM=  4.*X*SAX*(1.+2.*XX)/((QA-Q8)*A) 

QPP=  2.*SAX*QA/((QA-QB)**2) 

QPS=  SAX/ (QA-QB) 

QSP=-SAX/(QA-QB) 

QSS=  -2.*SAX*QB/ ( (QA-QB )**2) 

QPHI=  2.*(QN1/SAX  + QN2/SX) 

RETURN 

END 


GR009900 

GR009910 

GR009920 

GR009930 

GR009940 

GR009950 

GR009960 

GR009970 

GR009980 

GR009990 

GR010000 

GR010010 

GR010020 

GR010030 

GR010040 

GR010050 

GR010060 

GR010070 

GR010080 

GR010090 

GR010100 

GR0101 10 

GRO 101 20 

GR010130 

GR010140 

GR010150 

GR010160 

GR010170 

GR010180 

GR010190 

GR010200 

GR010210 

GR010220 

GR010230 

GR010240 

GR010250 

GR010260 

GR010270 

GR010280 

GR010290 

GR010300 

GR010310 

GR010320 

GR010330 

GR010340 

GR010350 

GR010360 

GR010370 

GR010380 

GR010390 

GR010400 

GR010410 

GR010420 

GR010430 

GRO 10440 


Q ***************************************************************** 

c 

INTEGER  FUNCTION  BINO(J1,J2) 

C 

Q ***************************************************************** 

c 

M = J1-2 

N = J2-2 

IF  ((M.EQ.N).OR.(N.EQ.O))  THEN 
BINO=1 
RETURN 

ELSE  IF  ((M.LT.O).OR.(N.LT.O).OR.CM.LT.N))  THEN 
BINO  = 0 
RETURN 

ELSE 
K=M-N 
IN  = 1 


IP  = M 
IQ  = 1 
DO  20  1=1, K 

IN  = IN* IP/ IQ 
IP  = IP-1 
IQ  = IQ  +1 
20  CONTINUE 

BINO  = IN 
RETURN 
ENDIF 
END 


► ★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★■AT**** 


C 


FUNCTION  DXM(D,M) 
C 


C 

C THE  FUNCTION  DXM(D,M)  IS 

C THE  SAME  AS  D**M.  BOTH  D AND  DXM  MUST  BE  DECLARED  IN 
C THE  MAIN  PROGRAM  AS  DOUBLE  PRECISION  VARIABLES. 

C 

DOUBLE  PRECISION  DXM,D 
IF  (M.EQ.O)  DXM=1.0D0 
IF  (M.GT.O)  GO  TO  11 
IF  (M.LT.O)  GO  TO  22 
RETURN 
C 

11  DXM=1.0DO 

DO  10  1=1, M 
DXM=DXM*D 
10  CONTINUE 

RETURN 
C 

22  DXM=1.0D0 

N=-M 

DO  20  1 = 1, N 
DXM=DXM/D 


GR010450 
GR010460 
GRO 10470 
GR010480 
GR010490 
GR010500 
GR010510 
GR010520 
GR010530 
GR010540 
GR010550 
GR010560 
GR010570 
GR010580 
GR010590 
GR010600 
GR010610 
GR010620 
GR010630 
GR010640 
GR010650 
GRO 10660 
GRO 10670 
GR0106S0 
GRO 10690 
GR010700 
GR010710 
GR010720 
GR010730 
GR010740 
GR010750 
GR010760 
GR010770 
GR010780 
GRO 10790 
GR010800 
GR010810 
GR010820 
GR010830 
GR010840 
GR010850 
GRO 10860 
GR010870 
GR010880 
GR010890 
GR010900 
GR010910 
GR010920 
GR010930 
GR010940 
GR010950 
GRO 10960 
GRO 10970 
GRO 10980 
GR010990 


20 

CONTINUE 

GR011000 

RETURN 

GR011010 

C 

GR011020 

END 

GR011030 

GR011040 

C 

GR011060 

SUBROUTINE  SFWAVE(ALPHA,XR, TTIME, INDEX, DISPL) 

GRO 11060 

C 

GR011070 

c** 

♦ ♦♦♦♦♦♦♦♦♦♦★★♦★★★★★★★★★★★★★★★★★★★★•Jrt***************************** 

GR011080 

c 

SFWAVE 

GR011090 

c 

SUBROUTINE  TO  COMPUTE  THE  SURFACE  DISPLACEMENT 

GRO 111  00 

c 

INPUTS: 

GR011110 

c 

ALPHA=RATIO  OF  SHEAR  WAVE  SPEED  TO  LONG.  WAVE 

SPEED 

GR011120 

c 

XR  =D I STANCE  BETWEEN  SOURCE  AND  DETECTOR 

GR011130 

c 

=ACTUAL  DIST./H 

GR011140 

c 

TTIME  =NOND IMENSIONALIZED  TTIME 

GR011150 

c 

=ACTUAL  TIME*SHEAR  WAVE  SPEED  / H 

GR011160 

c 

INDEX=SU8SCRIPT  OF  THE  GREEN'S  FUNCTION 

GR011170 

c 

e.g.  11,22,33,12,13  or  131,113,  etc. 

GR011180 

c 

GR011190 

c 

OUTPUT : 

GR011200 

c 

D I SPL=NOND IMENSIONALIZED  DISPLACEMENT 

GR011210 

c 

=ACTUAL  DISPL.*  PI  * SHEAR  MODULUS  * 

H / FORCE 

GR011220 

c 

NOTE: 

GR011230 

c 

THE  PARAMETERS  IN  COMMON  BLKS  0,1,  & 5 SHOULD 

BE  ENTERED 

GR011240 

c 

BEFORE  THE  FIRST  CALL  OF  THIS  SUBROUTINE, 

GR011250 

c 

I.  E.  CALL  INIT(ALPHA,XR, ZD, INDEX)  FIRST. 

GR011260 

c 

GR011270 

c 

SUBROUITNES  REQUIRED: 

GR011280 

c 

GR011290 

c 

DSFINT,DQR,QU,COMMON,RAYRT  & INTEGRATION  ROUTINES 

GR011300 

c 

GR011310 

c 

REMARK:  THE  COMPUTATION  IS  DONE  IN  Y-PLANE,  4tn  QUARD. 

ON  A 

GR011320 

c 

THREE  LEGGED  U-PATH  TO  REPLACE  THE  INTEGRATION 

FROM 

GR011330 

c 

-V  TO  V. 

GR011340 

c 

GR011350 

c 

GR011360 

c 

GR011370 

IMPLICIT  DOUBLE  PRECISION  (D) 

GR011380 

IMPLICIT  DOUBLE  COMPLEX  (Y) 

GR011390 

COMMON  /BLKO/DXD,DT,DZ,DA 

GR011400 

COMMON  /BLK1/KASE,M,L,K,N1,N2,D,DTHETA 

GR011410 

COMMON  /BLK5/YR, Y1 , Y2 

GR011420 

c 

GR011430 

N1=0 

GR011440 

N2=0 

GR011450 

DYR=OBLE(YR) 

GR011460 

DT=DBLE(TTIME) 

GR011470 

DZERO=O.ODO 

GR011480 

DMONEH=-0.5D0 

GR011490 

DSQRYR=DSQRT (DYR) 

GR011500 

RWSPD=DSQRYR 

GR011510 

c 

GR011520 

DPF=1.2D0/DBLE(ALPHA) 

GR011530 

DV=DXD/DT 

GR011540 

IF  (TTIME .LE . (ALPHA*XR) ) THEN 
DISPL=O.ODO 
RETURN 

ELSE 

DISPL=O.DO 

D I SPL=D I SPL+DSF I NG(DZERO, DMONEH , 1 ) 

2 +DSFING(DV,DPF,2) 

3 +DSFING(DMONEH,DZERO,3) 

END  IF 

END 

c 

DOUBLE  PRECISION  FUNCTION  DSFINGCDINIT,DEND, I PATH) 

C 

£************1Ht****************************************************** 

C 

IMPLICIT  DOUBLE  PRECISION  (D) 

C 

DOUBLE  PRECISION  EPS,R, E ,U(50,6) , FMIN, FMAX, F 
INTEGER  NINT,NMAX,KF, I FLAG 
LOGICAL  RST 
C 
C 

EXTERNAL  DQR1 ,DQR2,DQR3 
C 

EPS  = 1.D-8 
RST  = .FALSE. 

NMAX  =50 
NINT  = 1 
W(1,1)  = DINIT 
W(2,1)  = DEND 
IF  ( IPATH.EQ. 1 ) THEN 

CALL  DGLQ1 (DQR1,DINIT , DEND, EPS, R,E, NINT, RST, W,NMAX, FMIN, FMAX, 
+ KF, I FLAG) 

ELSE  IF  (IPATH.EQ. 2)THEN 

CALL  DGLQ1(DQR2,DINIT,DEND,EPS,R,E, NINT, RST,W,NMAX, FMIN, FMAX, 
+ KF, I FLAG) 

ELSE  IF  (IPATH.EQ. 3)  THEN 

CALL  DGLQ1 (DQR3, DINIT, DEND, EPS, R,E, NINT, RST, W, NMAX, FMIN, FMAX, 
+ KF, I FLAG) 

ELSE 

PAUSE  'DSFING:  IPATH  .NE.  1,2,  OR  3' 

END  IF 
C 

D2PEI=8.D0*DASIN(1 .DO) 

DSFING=-R/D2PEI 

RETURN 

END 

C******************************************************************* 

c 

DOUBLE  PRECISION  FUNCTION  DQRI(DL) 

C 

C****  *******  **************  ****************************************** 

IMPLICIT  DOUBLE  PRECISION  (D) 

EXTERNAL  DQR 


GR011550 

GR011560 

GR011570 

GR011580 

GR011590 

GR011600 

GR011610 

GRO 11620 

GR011630 

GR011640 

GR011650 

GR011660 

GR011670 

GRO 11680 

GR011690 

GR011700 

GR011710 

GR011720 

GR011730 

GRO 11740 

GR011750 

GR011760 

GR011770 

GRO 11780 

GRO 11790 

GR011800 

GR011810 

GR011820 

GRO 11830 

GR011840 

GR011850 

GR011860 

GR011870 

GRO 11880 

GR011890 

GR011900 

GRO 1 1910 

GR011920 

GR011930 

GRO 11940 

GR011950 

GR011960 

GR011970 

GRO 11980 

GRO 11990 

GR012000 

GR012010 

GR012020 

GR012030 

GR012040 

GR012050 

GR012060 

GR012070 

GR012080 

GR012090 


DQR1 =DQR(DL, 1 ) 

RETURN 

END 

(2  ******************************************************************* 

C 

DOUBLE  PRECISION  FUNCTION  DQR2(DL) 

C 

£******************************************************************* 


IMPLICIT  DOUBLE  PRECISION  (D) 
EXTERNAL  DQR 
DQR2=DQR(DL, 2) 

RETURN 

END 


C 

DOUBLE  PRECISION  FUNCTION  DQR3CDL) 
C 


IMPLICIT  DOUBLE  PRECISION  (D) 

EXTERNAL  DQR 
DQR3=DQR(DL,3) 

•RETURN 

END 

£★★★★★★★★★★★*★★*★★★★*★★*★**★★★*★★★★★*★★*★★*★★★★*★**★★***★**★★*★****★ 


C 


DOUBLE  PRECISION  FUNCTION  DQR(DL, IPATH) 
C 


C 


IMPLICIT  DOUBLE  COMPLEX  (Y) 

IMPLICIT  DOUBLE  PRECISION  (D) 

DIMENSION  YQD(3,3) , YUP(3,3) ,YQK(3) , YFEE(3) 

COMMON  /BLKO/DXD,DT,DZ,DA 

COMMON  /BLKl/KASE,M,L,KfN1,N2,D,DTHETA 

COMMON  /BLK2/YD, Y, YP, YQ, YR0,YRP, YRM, YDELTA,YPSQP1 , YQSQP1 , 

1 YPSQM1 , YQSQM1 , YSQRP1 , YSQRQ1 , YY, YA, YAA, YETA1 , YETA2 

2 , YDPHI 

COMMON  /BLK3/YI , YONE, YTWO, YTHREE, YFOUR, YEIGHT, YZERO 
C 

c 

DV=DXD/DT 

IF  (IPATH. EQ.1)  THEN 
Y=  DCMPLX(  DV,DL) 

ELSE  IF  (IPATH. EQ. 2)  THEN 
Y=  DCMPLX(DL, -0.5D0) 

ELSE  IF  (IPATH. EQ. 3)  THEN 
DPF=1 .2D0/DA 
Y=  DCMPLX(DPF,DL) 

ELSE 

PAUSE  'DRQ:  IPATH  ?' 

END  IF 

CALL  SFCOM(DL) 

CALL  QD(YQD, IDUMMY) 

CALL  QU(YUP, IDUMMY) 

C 


GR012100 

GR012110 

GR012120 

GR012130 

GR012140 

GR012150 

GR012160 

GR012170 

GR012180 

GR012190 

GR012200 

GR012210 

GR012220 

GR012230 

GR012240 

GR012250 

GR012260 

GR012270 

GR012280 

GR012290 

GR012300 

GR012310 

GR012320 

GR012330 

GR012340 

GR012350 

GR012360 

GRO 12370 

GR012380 

GR012390 

GR012400 

GR012410 

GR012420 

GR012430 

GR012440 

GR012450 

GRO 12460 

GR012470 

GR012480 

GR012490 

GR012500 

GRO 1 25 1 0 

GR012520 

GR012530 

GR012540 

GR012550 

GR012560 

GR012570 

GR012580 

GR012590 

GR012600 

GR012610 

GR012620 

GR012630 

GR012640 


YFEEC 1 )=YAA/ (YP*YD) 

GR012650 

YFEE(2)=Y0NE/(YQ*YD) 

GRO 12660 

YFEE(3)=Y0NE/(YQ*YD) 

GR012670 

c 

GRO 12680 

YQR  = YZERO 

GR012690 

c 

GR012700 

IF  (K.EQ.O)  THEN 

GR012710 

YQK(1 )=YONE 

GR012720 

YQK(2)=YONE 

GR012730 

YQK(3)=YONE 

GRO 12740 

ELSE  IF  (K.EQ.1)  THEN 

GR012750 

YQK(1 )=YETA1/Y 

GRO 12760 

YQK(2)=YETA1/Y 

GR012770 

YQK(3)=YETA1/Y 

GR012780 

ELSE  IF  OC.EQ.2)  THEN 

GRO 12790 

YQK(  1 )=YETA2/Y 

GR012800 

YQK(2)=YETA2/Y 

GR012810 

YQK(3)=YETA2/Y 

GR012820 

ELSE  IF  (K.EQ.3)  THEN 

GR012830 

YQK( 1 )=- YP/Y 

GRO 12840 

YQK(2)=- YQ/Y 

GR012850 

YQK(3)=- YQ/Y 

GR012860 

ELSE 

GRO 12870 

PAUSE  ' FUUCTION  DOR:  WRONG  INDEX,  < .NE.  0,1,2,  OR  3' 

GRO 12880 

END  IF 

GR012890 

DO  100  1=1,3 

GR012900 

YQR=YQR+YQD(M, I )*YUP(L, I )*YQK( I ) *Y  FEE  C I ) 

GR012910 

100  CONTINUE 

GRO 12920 

c 

GR012930 

IF  ( ( I PAT H.EQ.1).OR.(I PATH . EQ . 3 ) ) THEN 

GR012940 

YQR=YQR*YI/CDSQRT (DCMPLX(DV*DV,0.D0) - YY) 

GR012950 

ELSE 

GRO 12960 

YQR=YQR/CDSQRT (DCMPLX(DV*DV, 0.D0) - YY) 

GRO 12970 

END  IF 

GRO 12980 

DQR=DBLE(YQR) 

GR012990 

c 

GR013000 

RETURN 

GR013010 

END 

GR013020 

^★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★*****************************Hr****** 

GR013030 

c 

GR013040 

SUBROUTINE  SFCOM(DL) 

GR013050 

c 

GR013060 

0***r***********************1***************************************** 

GR013070 

c 

GR013080 

c 

TO  COMPUTE  PARAMETERS  DEPENDING  ON  Y AND  TO  PUT  IN 

GR013090 

c 

A COMMON  BLK2 

GR013100 

c 

GR0131 10 

IMPLICIT  DOUBLE  COMPLEX  (Y) 

GR013120 

IMPLICIT  DOUBLE  PRECISION  (D) 

GR013130 

DOUBLE  PRECISION  ETA1 ,ETA2SQ,DD,AA 

GR013140 

COMMON  /BLK0/DXD,DT,DZ,DA 

GR013150 

COMMON  /BLK1/KASE,M,L,K,N1 ,N2,D,DTHETA 

GR013160 

COMMON  /BLK2/YD, Y, YP, YQ, YRO, YRP, YRM, YDELTA, YPSQP1 , YQSQP1 , 

GR013170 

1 YPSQM1 , YQSQM1 , YSQRP1 , YSQRQ1 , YY , YA , YAA,  YETA1 , YETA2 

GR013180 

2 , YDPHI 

GR013190 

COMMON  /BLK3/YI , YONE, YTWO, YTHREE, Y FOUR, YE I GHT, YZERO 
C 

AA=DA*DA 

DV=DXD/DT 

DN1=0BLE(REAL(N1 )) 

DN2=DBLE(REAL(N2)> 

YZERO=OCMPLX( 0.00,0. DO) 

YI=DCMPLX(0.D0, 1.0D0) 

YAA=DCMPLX(AA, 0.D0) 

YA=DCMPLX(DA,0.D0) 

YONE=DCMPLX(1. 000,0.000) 

YTWO=DCMPLX(2.DO,O.DO) 

YFOUR=DCMPLX(4. 000,0. DO) 

C 

YD=- Y*DCMPLX(DT,0.D0) 

YETA1  = -Y/OCMPLX(DV,O.DO) 

YETA2=CDSQRT(YONE-YETA1*YETA1 ) 

IF  (OBLE(Y) .LT.O.DO)  YETA2=-YETA2 
YY=Y*Y 

YP=CDSQRT ( YAA*YY - YONE ) 

YQ=CDSQRT<YY-YONE) 

YPSQP1=YAA*YY 
YQSQP1=YY 
YQSQM1=YQ*YQ- YONE 
YPSQM1=YP*YP- YONE 
C 

YSQRP1=Y*YA 

YSQRQ1=Y 

YDELTA=YQSQM1*YQSQM1+YF0UR*YP*YQ 
YM=(YQSQM1*YQSQM1 - YFOUR*YP*YQ) 

YRO=YM/YDELTA 

YRP=* (YF0UR*YP*YQSQM1 )/(YA*YDELTA) 

YRM=(YRO*YRO-YONE)/YRP 

C 

RETURN 

END 

C******************************************************************* 


SUBROUTINE  QU(YUP, IDUMMY) 


£★★*★*★*★★**★★*★★★**★**★★**★****★*★*★*★★******★★**★★★*★★★*★★*★★★*★★★ 


C 

C TO  COMPUTE  QU(I,J)  BASED  ON  EQ.  6.4 
C 

IMPLICIT  DOUBLE  COMPLEX  (Y) 

COMMON  /BLKO/DXD , DT , DZ , DA 

COMMON  /BLK1/KASE,M,L,K,N1,N2,D,DTHETA 

COMMON  /BLK2/YD,Y, YP,YQ,YR0, YRP, YRM, YDELTA, YPSQP1 , YQSQP1 , 

1 YPSQM1 , YQSQM1 , YSQRP1 , YSQRQ1 , YY , YA , YAA , YETA1 , YETA2 

2 , YDPHI 

COMMON  /BLK3/Y I , YONE , YTWO , YTHREE , YFOUR , YE I GHT , YZERO 
C 

DIMENSION  YUP(3,3) 

C 

IDUMMY=0 


GR013200 
GR013210 
GR013220 
GRO 13230 
GR013240 
GR013250 
GR013260 
GR013270 
GR013280 
GR013290 
GR013300 
GR013310 
GR013320 
GR013330 
GR013340 
GR013350 
GR013360 
GR013370 
GR013330 
GR013390 
GR013400 
GR013410 
GR013420 
GR013430 
GR013440 
GR013450 
GR013460 
GRO 13470 
GR013480 
GR013490 
GR013500 
GR013510 
GR013520 
GR013530 
GR013540 
GR013550 
GR013560 
GR013570 
GR013580 
GR013590 
GR013600 
GR013610 
GR013620 
GR013630 
GR013640 
GR013650 
GR013660 
GR013670 
GRO 13680 
GR013690 
GR013700 
GR013710 
GR013720 
GR013730 
GR013740 


YUP(1 , 1 )=YETA1/YSQRP1 
YUP(1,2)=-YETA1*YQ/YSQRQ1 
YUP< 1 , 3)=- YETA2 
YUP(2, 1 )=YETA2/YSQRP1 
YUP(2,2)=- YETA2*YQ/YSQRQ1 
YUP(2,3)=YETA1 
YUP(3f 1 )=-YP/YSQRP1 
YUP(3f 2)=- Y0NE/YSQRQ1 
YUP(3,3)=YZERO 
C 

RETURN 

END 

C******************************************************************* 

C 

SUBROUTINE  RAYRT(YR,YI 1 ,YI2, ALPHA) 

C 

Q******************************************************************* 

c 

C TO  SOLVE  RAYLEIGH  EQUITION 
C 

C INPUT:  ALPHA  = SHEAR  WAVE  SPEED  / LONG.  WAVE  SPEED 
C 

C OUTPUT:  YR, YI 1 , YI2  ARE  THE  ROOTS  TO 
C 

C Y*Y*Y-8*Y*Y+8*(3-2*AA)*Y- 16*( 1 - AA)  = 0 
C 

DOUBLE  COMPLEX  YR,YI1,YI2 

DOUBLE  PRECISION  A,AA,RT1 ,RT2,RT3 

DIMENSION  A(4),RT1(2),RT2(2),RT3(2) 

C 

I F( (ALPHA. LE.O.). OR. (ALPHA. GE. 1.0))  THEN 
WRITE(6,*) 'RAYRT : WRONG  VALUE  OF  ALPHA' 

WRITE(6,*)  ALPHA 

WRITE (6,*) 'RAYRT : COMPUTATION  PAUSED' 

PAUSE 
END  IF 
C 

AA=ALPHA*ALPHA 
TOL=  1.E-37 
A( 1 ) = 1.0 
A(2)  = -8.0DO 
A(3)  = 8.*(3. -2.*AA) 

A(4)  = -16.*(1 . -AA) 

C 

CALL  CUB  I C( A , RT 1 , RT2 , RT3 , TOL , NOROOT ) 

C 

C IF  (ALPHA.LT. 0.567)  THEN 

C 

C ONLY  ONE  REAL  ROOT 

C 

C ELSE 

C 

C ALL  THREE  ROOTS  ARE  REAL 

C 

IF  ( (RT1 ( 1 ) . LT . RT2( 1 ) ) .AND . (RT1 ( 1 ) . LT.RT3( 1 ) ) )THEN 


GR013750 
GRO 13760 
GR013770 
GR013780 
GRO 13790 
GR013800 
GR013810 
GR013820 
GR013830 
GR013840 
GR013850 
GR013860 
GRO 13870 
GR013880 
GR013890 
GR013900 
GR013910 
GR013920 
GR013930 
GR013940 
GR013950 
GR013960 
GR013970 
GR013980 
GR013990 
GR014000 
GR014010 
GR014020 
GR014030 
GR014040 
GR014050 
GR014060 
GR014070 
GR014080 
GR014090 
GR014100 
GR0141 10 
GRO 141 20 
GR014130 
GR014140 
GR014150 
GR014160 
GR014170 
GR014180 
GR014190 
GR014200 
GR014210 
GR014220 
GR014230 
GR014240 
GR014250 
GR014260 
GR014270 
GR014280 
GR014290 


YR=DCMPLX(RT1 (1),RT1(2)) 

YI1=DCMPLX(RT2(1 ),RT2(2)) 

YI2=DCMPLX(RT3(1),RT3(2)) 

ELSEIF  ((RT2(1).LT.RT1(1)) .AND . (RT2( 1 ) . LT .RT3( 1 ) ) ) THEN 
YI1=DCMPLX(RT1(1),RT1(2)) 

YR  =DCMPLX(RT2(1 ),RT2(2)) 

YI2=DCMPLX(RT3(1 ),RT3(2)) 

ELSEIF  ( (RT3( l).LT.RTI(D) .AND . (RT3(  1 ) . LT .RT2(  1 ) ) )THEN 
YI2=DCMPLX(RT1 ( 1),RT1(2)) 

YI 1=0CMPLX(RT2( 1 ) ,RT2(2) ) 

YR  =DCMPLX(RT3( 1 ) , RT3(2) ) 

END  IF 
RETURN 
END 


C 

SUBROUT  I NE  T I MEARR I ( A , R , ZD , NRAY , RCN , TA , CN ) 

C 

c******************************************************************* 

c 

C 9/27/82 

C 

C 

DIMENSION  Z1(4),Z2(4),  TA(NRAY),CN(NRAY,3),RCN(NRAY,3) 

DOUBLE  COMPLEX  YR,YI1,YI2 
I TER= 1 00 
NRAY=NRAY 
1=0 
LM=1 

C Z = ZD+0.5  ; ZD  refers  to  origin  at  the  center  of  the  plate 
Z = ZD+0.5 
Z1(1)=0.0 
Z2(1 )=Z 
Z1 (2)=Z 
Z2(2)=0.0 
Z1(3)=-Z 
Z2(3)=0.0 
Z1(4)=0.0 
Z2(4)=-Z 
C 

IF  (Z.EQ.1.0)  Z2(1)=0.0 
OREVEN=FLOAT(INT(Z))*0.5 
DO  100  J=1 ,21 
DO  200  K=1 ,12 
JJ=J-1 
KK=K- 1 

IF  ( (K.EQ.I).AND.(J.EQ.I)  ) THEN 
IF(Z.EQ.I.O)  THEN 

C THREE  SURFACE  WAVE  ARRIVALS: 

1 = 1+1 

CN( I , 1 )=- 1 . 

CN( I ,2)=0. 

CN( I ,3)=R*A 
TA( I )=CN( I ,3) 

1 = 1+1 


GR014300 
GR014310 
GR014320 
GR014330 
GR014340 
GR014350 
GR014360 
GR014370 
GR014380 
GRO 14390 
GR014400 
GR014410 
GR014420 
GRO 14430 
GR014440 
GR014450 
GR014460 
GR014470 
GR014480 
GR014490 
GR014500 
GR014510 
GR014520 
GR014530 
GR014540 
GR014550 
GR014560 
GR014570 
GR014580 
GR014590 
GR014600 
GR014610 
GR014620 
GR014630 
GR014640 
GR014650 
GR014660 
GR014670 
GR014680 
GR014690 
GR014700 
GR014710 
GR014720 
GR014730 
GR014740 
GR014750 
GRO 14760 
GR014770 
GR014730 
GR014790 
GR014800 
GR014810 
GR014820 
GR014830 
GR014840 


CN( I , 1 )=- 1 . 

CN(I,2)=-1. 

CN( I ,3)=R 
TA( I )=CN( 1,3) 

CALL  RAYRT (YR,YI1,YI2,A) 

1 = 1 + 1 

CN( I , 1 )=0. 

CN( I ,2)=0. 

CN  < 1 ,3)=R/DSQRT(DBLE(YR)) 

TA( I )=CN( I ,3) 

END  IF 
GOTO  200 
ENDIF 

IF  ( (Z.EQ.O.O).OR.(Z.EQ.I.O)  ) LM=4 
IF  ( ((FLOAT( JJ+KK))/2.0-OREVEN) .EQ. FLOAT ( ( J J+KK)/2)  ) 
1 GO  TO  200 

DO  300  LL=1,4,LM 
C1=JJ+Z1(LL) 

C2=KK+Z2(LL) 

IF  ( (C1.LT.0.0)  .OR.  (C2.LT. 0.0)  ) GO  TO  300 
1 = 1+1 

IF  (I.GT.NRAY)  GO  TO  99 

CALL  TARRV  (A,C1 ,C2,R, I TER , ERROR, T) 

IF  (ERROR. GT.1.E-6)  WRITE  (6,603)  I ,C1 ,C2, ERROR 
CN( I , 1 )=C1 
CN( I ,2)=C2 
CN( I ,3)=T 
TA( I )=T 

IF  (CI.GT.O.)GO  TO  300 
APJ=A*C2/SQRT ( 1 . - A*A) 

IF(R.LE.APJ)  GO  TO  300 
C2=-C2 

CALL  TARRV  (A,C1 ,C2,R, I TER, ERROR, T) 

1 = 1+1 

IF  (I.GT.NRAY)  GO  TO  99 
CN(  1 , 1 )=C1 
CN( I ,2)=C2 
CN( I ,3)=T 
TA( I )=T 

300  CONTINUE 

200  CONTINUE 
100  CONTINUE 
C 


99  CONTINUE 

CALL  RSORT(CN,TA,RCN,NRAY,3,0) 

C 

603  FORMAT(2X, 'ERROR  IS  LARGER  THAN  1.E-6  • , I3,3F10.5) 
RETURN 
END 


C 

SUBROUTINE  TARRV(A,C1 ,C2,R, ITER,ERR,TARR) 

C 

C******************************************************************* 


C 


GR014850 

GRO 14860 

GRO 14870 

GR014880 

GRO 14890 

GR014900 

GR014910 

GR014920 

GR014930 

GR014940 

GR014950 

GR014960 

GR014970 

GR014980 

GR014990 

GR015000 

GR015010 

GR015020 

GR015030 

GR015040 

GR015050 

GR015060 

GR015070 

GR015080 

GR015090 

GRO 1 5 1 00 

GR0151 10 

GR015120 

GR015130 

GRO 151 40 

GR015150 

GR015160 

GR015170 

GR015180 

GR015190 

GR015200 

GR015210 

GR015220 

GR015230 

GR015240 

GR015250 

GR015260 

GR015270 

GR015280 

GR015290 

GR015300 

GR015310 

GR015320 

GR015330 

GR015340 

GR015350 

GR015360 

GR015370 

GRO 15380 

GR015390 


IMPLICIT  DOUBLE  PRECISION  (D) 

GR015400 

c 

GR015410 

NUM  = 0 

GR015420 

DTI =0.1570796326794901 

GR015430 

DT2=0.0D0 

GR015440 

DA=D8LE(A) 

GR015450 

DAA=DA*DA 

GR015460 

DN1=DBLE(C1 ) 

GR015470 

DN2=DBLE(C2) 

GR015480 

DR=DBLE(R) 

GR015490 

IF  (C1.EQ.0.0)  GO  TO  2 

GR015500 

IF  (C2.EQ.0.0)  GO  TO  3 

GR015510 

c 

GR015520 

1 

CONTINUE 

GR015530 

NUM  = NUM  + 1 

GR015540 

DT=(DT1+DT2)/2.0D0 

GR015550 

DAASS=DAA*DSIN(DT)*DSIN(DT) 

GR015560 

DAASS=1 .ODO-DAASS 

GR015570 

DCS2=DSQRT (DAASS) 

GR015580 

DRT  = DN1*DTAN(DT)  + DN2*DA*DSIN(DT)/DCS2 

GR015590 

c 

URITE  (6,20)  NUM,DT,DRT,C1 ,C2 

GR015600 

IF  (NUM. GE. ITER)  GO  TO  10 

GR015610 

IF  (DABS(DRT-DR) .LE.5.D- 11)  GO  TO  10 

GR015620 

IF  (DRT.LE.DR)  DT2=DT 

GR015630 

IF  (DRT.GT.DR)  DT1=DT 

GR015640 

GO  TO  1 

GR015650 

c 

GR015660 

c 

20 

FORMAT  ('NUM=  ■ , 13, ' T=  * ,D17. 10, 1 RT=  ' ,017. 10,2(2X, F4. 0) ) 

GR015670 

10 

ITER=NUM 

GR015680 

DTARR=DA*DN1/DCOS(DT)  + DN2/DCS2 

GR015690 

DERR=DABS(DR-DRT) 

GR015700 

ERR=DERR 

GR015710 

TARR=DTARR 

GR015720 

RETURN 

GR015730 

c 

GR015740 

2 

IF  (C2.LT.0.0)  GO  TO  4 

GR015750 

DTARR=DSQRT (DR*DR  + DN2*DN2) 

GR015760 

TARR=DTARR 

GR015770 

ERR=0.0 

GR015780 

RETURN 

GR015790 

c 

GR015800 

3 

DTARR=DSQRT (DR*DR+DN1*DN1 )*DA 

GR015810 

TARR=DTARR 

GR015820 

ERR=0.0 

GR015830 

RETURN 

GR015840 

c 

GR015850 

4 

DTARR=DR*DA  - DN2*DSQRT ( 1 .DO-DAA) 

GR015860 

c 

NOTE  HERE  DN2  < 0 

GR015870 

TARR=DTARR 

GR015880 

ERR=0.0 

GR015890 

RETURN 

GR015900 

c 

GR015910 

END 

GR015920 

Q******************************************************************* 

GR015930 

c 

GR015940 

SUBROUT I NE  CUB I C( A , RT1 , RT2 , RT3 , TOL , NOROOT ) 

C 

£★**★★★★★★★*★***★★★★★★***★**★★*★★***★★★***★★**★***★**★*****★*★*★★★** 

C 

C THIS  SUBROUTINE  FINDS  THE  ROOTS  OF  A CUBIC  EQUITION 
C 

C INPUTS: 

C A COEFFICIENTS,  A(4) 

C 

C OUTPUTS: 

C RT1 ,RT2,RT3  THREE  ROOTS 

C RT1 (1 )=REAL  PART  OF  RT1 , ETC. 

C 

C THE  ROUTINE  CALLS  LINCNG,AND  QUADRT 

C 

DOUBLE  PRECISION  A,B,C,RT1 ,RT2,RT3,ZEROfX,Y,ZZ,SQT3,RC27fANG 
DIMENSION  A(4) ,RT1 (2) ,RT2(2) , RT3(2) , B(2) ,C(4) 

C 

ZERO=TOL/10. 

SQT3=DSQRT (3.0DO) 

RC27=1. 0/27.0 
C 

IF  (DABS(A( 1 ) ) -ZERO)  7,7,12 
7 CALL  QUADRT  (A(2),RT1 ,RT2, TOL, NOROOT) 

RETURN 

C 

12  NOROOT =3 
IT=0 
C 

DO  1 1=2,4 

A( I )=A( I )/A( 1 ) 

1  CONTINUE 
A( 1 )=1 .0 
C 

IF  (DABS(A(2) ) .LE .ZERO)  GOTO  2 
NDA=3 
B ( 1 )=1 .0 

B ( 2 )= - (A(2)/3.0) 

CALL  LINCNG(A,NDA,B,C) 

I T=  1 
GO  TO  4 

2 DO  3 1=1,4 

C( I )=A( I ) 

3 CONTINUE 
C 

4 X = C(4)*C(4)*0.25+(C(3)**3)*RC27 

IF  (X.LT.0.0)  GO  TO  9 
X = DSQRT(X) 

Y = -(C(4)*0.5) 

1=1 

RTI(I)  =Y+X 
C 

5 N = 0 

I F C RT 1 ( I ) .LT.0.0)  N=1 
IF  (DABS(RTI(D).LE.ZERO)  GO  TO  6 


GR015950 

GR015960 

GR015970 

GR015980 

GR015990 

GR016000 

GR016010 

GR016020 

GR016030 

GR016040 

GR016050 

GRO 16060 

GR016070 

GRO 16080 

GR016090 

GR016100 

GR016110 

GR016120 

GR016130 

GR016140 

GR016150 

GR016160 

GR016170 

GR016180 

GR016190 

GR016200 

GR016210 

GR016220 

GR016230 

GR016240 

GR016250 

GR016260 

GR016270 

GR016280 

GR016290 

GR016300 

GR016310 

GR016320 

GR016330 

GR016340 

GR016350 

GRO 16360 

GR016370 

GR016380 

GR016390 

GR016400 

GR016410 

GR016420 

GR016430 

GR016440 

GR016450 

GRO 16460 

GR016470 

GR016480 

GR016490 


RT 1 ( I ) = DABS(RT1(I))**(1./3.) 

GR016500 

IF  (N.EQ.1)  RT1(I)  = -RTI(I) 

GR016510 

6 

IF  (I.EQ.2)  GOTO  8 

GR016520 

I = 2 

GR016530 

RT 1 C I ) = Y-X 

GR016540 

GO  TO  5 

GR016550 

c 

GR016560 

8 

RT2(2)  = ( (RT1 ( 1 ) - RT1(2))*0.5)*SQT3 

GR016570 

RT1 ( 1 ) = RT1 (1 ) + RT1 (2) 

GR016580 

RT 1 < 2 > = 0.0 

GR016590 

RT2(1 ) = - RT 1 ( 1 )*0.5 

GR016600 

RT3( 1 ) = RT2(1 ) 

GR016610 

RT3C2)  = -RT2(2) 

GR016620 

GO  TO  11 

GR016630 

c 

GRO 16640 

9 

ZZ  = DABS(C(3>) 

GR016650 

X = -(C(4)*0.5)/DSQRT((ZZ**3)*RC27  ) 

GR016660 

ANG  = DACOS(X) 

GR016670 

Y = 2.0*(DSQRT (ZZ/3.0) ) 

GR016680 

ANG=ANG/3.0 

GRO 16690 

RT1 ( 1 ) = Y*DCOS(ANG) 

GR016700 

RT2C1)  = Y*DCOS(ANG  + 2. 094395 1024D0) 

GR016710 

RT3(1 ) = Y*DCOS(ANG  + 4. 1887902047D0) 

GR016720 

RT1 (2)  = 0.0 

GR016730 

RT2(2)  = 0.0 

GR016740 

RT3(2)  = 0.0 

GR016750 

c 

GRO 16760 

11 

IF  (IT.EQ.1)  THEN 

GRO 16770 

RT1 ( 1 ) = RT 1(1)  + B(2) 

GRO 16780 

RT2( 1 ) = RT2(1 ) + B(2) 

GR016790 

RT3(1 ) = RT3( 1 ) + B(2) 

GR016800 

ELSE 

GR016810 

END  IF 

GR016820 

RETURN 

GR016830 

END 

GR016840 

£★★***•**★**★*★**★ *★*★*★★★*★**★★***★*★★***★★****★★★★★**★★★★*★★★★*★★★★ 

GR016850 

c 

GR016860 

SUBROUTINE  LINCNGCA, NDA,B,C) 

GR016870 

c 

GR016880 

£****★************************************************************** 

GR016890 

c 

GR016900 

c 

TO  MAKE  A LINEAR  CHANGE  OF  VARIABLES 

GR016910 

c 

IN  A GIVEN  POLYNOMIAL 

GR016920 

c 

GR016930 

DOUBLE  PRECISION  A,B,C,BIPWR 

GR016940 

DIMENSION  A(1 ), B(2), C(1 ) 

GR016950 

c 

GR016960 

BIPWR  = 1.0D0 

GR016970 

NZ  = NDA  +1 

GR016980 

c 

GR016990 

DO  4 1=1 ,NZ 

GR017000 

C(I)  = A(I) 

GR017010 

4 

CONTINUE 

GR017020 

c 

GR017030 

6 

DO  8 1=2, NZ 

GR017040 

C(I)  = C(I)  + CC 1-1 )*B(2) 

GR017050 

CONTINUE 

GROI 7060 

GROI 7070 

C(NZ)  = C(NZ)*BIPWR 

GROI 7080 

NZ  = NZ  - 1 

GROI 7090 

BIPWR  = BIPWR*B( 1 ) 

GR017100 

IF  (NZ.GT.1)  GO  TO  6 

GR0171 10 

C(1)  = C( 1 )*BIPWR 

GR017120 

RETURN  GR017130 

END  GR017140 

£**■*■***  *■**■*  ■**-**********+***■***  ******  ********************************  GR01 71 50 

C GR017160 

SUBROUTINE  DGLQ1(F,A,B,EPSf RjE^INTjRST.W, NMAX,FMIN,FMAX,KF, I FLAG)GR017170 
C GR017180 

q** *************************************  ******* *★*★★*★★★★★★★★★★★*★★★  grOI 71 90 

C GROI 7200 

C GR017210 


C 

*** 

ALL  REAL  VARIABLES  ARE  TYPED  DOUBLE  PRECISION 

GR017220 

C 

★★★ 

9/20/82  NH 

GR017230 

C 

GROI 7240 

C 

GROI 7250 

C 

DGLQ1 

IS  A SUBROUTINE  FOR  THE  AUTOMATIC  EVALUATION 

GROI 7260 

C 

OF  DEFINITE  INTEGRALS  OF  A USER  DEFINED  FUNCTION 

GROI 7270 

C 

OF  ONE  VARIABLE  PROVIDING  FLEXIBLE  USAGE. 

GR017280 

C 

GROI 7290 

C 

FOR  AN  EASY  TO  USE  VERSION  SEE  SUBROUTINE  DGLQ. 

GR017300 

C 

GR017310 

C 

CAPABILITIES  OF  DGLQ1  (IN  ADDITION  TO  THOSE  OF  DGLQ) 

GR017320 

C 

INCLUDE: 

GR017330 

C 

ABILITY  TO  RESTART  A CALCULATION  TO  GREATER 

GR017340 

C 

ACCURACY  WITHOUT  PENALTY... 

GR017350 

c 

ABILITY  TO  SPECIFY  AN  INITIAL  PARTITION  OF 

GR017360 

c 

THE  INTEGRATION  INTERVAL... 

GROI 7370 

c 

ABILITY  TO  INCREASE  THE  WORK  SPACE  TO  HANDLE 

GROI 7380 

c 

MORE  DIFFICULT  PROBLEMS... 

GR017390 

c 

OUTPUT  OF  LARGEST/SMALLEST  INTEGRAND  VALUE  FOR 

GROI 7400 

c 

APPLICATIONS  SUCH  AS  GRAPHING... 

GR017410 

c 

GROI 7420 

c 

A R G 

UMENTS  IN  THE  CALL  SEQUENCE 

GROI 7430 

c 

GROI 7440 

c 

F 

(INPUT)  THE  NAME  OF  YOUR  INTEGRAND  FUNCTION. 

GR017450 

c 

THIS  NAME  MUST  APPEAR  IN  AN  EXTERNAL  STATEMENT 

GROI 7460 

c 

IN  ANY  PROGRAM  WHICH  CALLS  DGLQ1. 

GROI 7470 

c 

YOU  MUST  WRITE  F IN  THE  FORM 

GROI 7480 

c 

FUNCTION  F(X) 

GROI 7490 

c 

F=(EVALUATE  INTEGRAND  AT  THE  POINT  X) 

GROI 7500 

c 

RETURN 

GR017510 

c 

END 

GROI 7520 

c 

A 

GR017530 

c 

B 

(INPUT)  ENDPOINTS  OF  INTEGRATION  INTERVAL 

GR017540 

c 

EPS 

(INPUT)  ACCURACY  TO  WHICH  THE  INTEGRAL  IS  TO  BE  CALCULAGR017550 

c 

DGLQ1  WILL  TRY  TO  ACHIEVE  RELATIVE  ACCURACY 

, GR017560 

c 

E.G.  SET  EPS=.01  FOR  2 DIGITS,  .001  FOR  3, 

ETGR017570 

c 

R 

(OUTPUT)  THE  ESTIMATE  OF  THE  INTEGRAL 

GROI 7580 

c 

E 

(OUTPUT)  THE  ESTIMATE  OF  THE  ABSOLUTE  ERROR  IN  R. 

GR017590 

NIMT  (INPUT  GR017600 

OUTPUT)  GR017610 

AS  AN  OUTPUT  QUANTITY,  NINT  GIVES  THE  GR017620 

NUMBER  OF  SUBINTERVALS  IN  THE  FINAL  GR017630 

PARTITION  OF  [A,B] . GR017640 

AS  AN  INPUT  QUANTITY,  NINT  MUST  8E  SET  TO  GR017650 

THE  NUMBER  OF  SUBINTERVALS  IN  THE  INITIAL  GR017660 

PARTITION  OF  [A,B] . FOR  MOST  PROBLEMS  GR017670 

THIS  IS  JUST  1,  THE  INTERVAL  [A,B]  ITSELF.  GR017680 
NINT  IS  USEFUL  IF  YOU  WOULD  LIKE  TO  HELP  GRO 17690 

DGLQ1  LOCATE  A DIFFICULT  SPOT  ON  [A,B] . GR017700 

IN  THIS  REGARD  NINT  IS  USED  ALONG  GR017710 

WITH  THE  ARRAY  W (SEE  BELOW).  IF  YOU  SET  GR017720 

N I NT = 1 IT  IS  NOT  NECESSARY  TO  BE  CONCERNED  GR017730 

WITH  W,  EXCEPT  THAT  IT  MUST  BE  DIMENSIONED. . .GR01 7740 
AS  AN  EXAMPLE  OF  MORE  GENERAL  APPLICATIONS,  GR017750 
IF  [A,B] = [0, 1]  BUT  THE  INTEGRAND  JUMPS  AT  0.3GR017760 
IT  WOULD  BE  WISE  TO  SET  NINT=2  AND  THEN  SET  GR017770 

W(1,1)=0.0  (LEFT  ENDPOINT)  GR017780 

W(2, 1 )=0.3  (SINGULAR  POINT)  GR017790 

W(3, 1 )=1 .0  (RIGHT  ENDPOINT)  GR017800 

IF  YOU  SET  NINT  GREATER  THAN  1,  BE  SURE  TO  GR017810 
CHECK  THAT  YOU  HAVE  ALSO  SET  GR017820 

W( 1 , 1 )=A  AND  W(NINT+1 , 1 )=B  GR017830 

RST  (INPUT)  A LOGICAL  VARIABLE  (E.G.  TRUE  OR  FALSE)  GR017840 

SET  RST=. FALSE.  FOR  INITIAL  CALL  TO  DGLQ1  GR017850 

SET  RST=.TRUE.  FOR  A SUBSEQUENT  CALL,  GR017860 

E.G.  ONE  FOR  WHICH  MORE  ACCURACY  IS  GR017870 

DESIRED  (SMALLER  EPS).  A RESTART  ONLY  GR017880 

MAKES  SENSE  IF  THE  PRECEDING  CALL  RETURNED  GR017890 
WITH  A VALUE  OF  I FLAG  (SEE  BELOW)  LESS  THAN  GR017900 
ON  A RESTART  YOU  MAY  NOT  CHANGE  THE  VALUES  0GR017910 
ANY  ARGUMENTS  IN  THE  CALL  SEQUENCE,  EXCEPT  EGR017920 
W(NMAX,6)  W IS  AN  ARRAY  USED  FOR  SCRATCH  STORAGE  BY  DGLQ1 .GR017930 

YOU  MUST  INCLUDE  A DIMENSION  STATEMENT  IN  GR017940 

YOUR  CALLING  PROGRAM  TO  ALLOCATE  THIS  STORAGE.  GR017950 

THIS  SHOULD  BE  OF  THE  FORM  GR017960 

DIMENSION  W(NMAX,6)  GR017970 

WHERE  NMAX  IS  AN  INTEGER.  AN  ADEQUATE  VALUE  OF  GR017980 
NMAX  IS  50.  GR017990 

NMAX  (INPUT)  AN  INTEGER  EQUAL  TO  THE  FIRST  SUBSCRIPT  IN  THE  GR018000 
DIMENSION  STATEMENT  FOR  THE  ARRAY  W.  THIS  IS  GR018010 
ALSO  EQUAL  TO  THE  MAXIMUM  NUMBER  OF  SUB  INTERVAL GR0 18020 
PERMITTED  IN  THE  INTERNAL  PARTITION  OF  [A,B].  GR018030 

A VALUE  OF  50  IS  AMPLE  FOR  MOST  PROBLEMS.  GR018040 

FMIN  GR018050 

FMAX  (OUTPUT)  THE  SMALLEST  AND  LARGEST  VALUES  OF  THE  I NTEGRANGR018060 
WHICH  OCCURRED  DURING  THE  CALCULATION.  THE  GR018070 

ACTUAL  INTEGRAND  RANGE  ON  [A,B]  MAY,  OF  COURSGRO 18080 
BE  GREATER  BUT  PROBABLY  NOT  BY  MORE  THAN  10X.GR018090 
KF  (OUTPUT)  THE  ACTUAL  NUMBER  OF  INTEGRAND  EVALUATIONS  USE0GR018100 

BY  DGLQ1  TO  APPROXIMATE  THIS  INTEGRAL.  KF  GR018110 
WILL  ALWAYS  BE  AT  LEAST  30.  GR018120 

I FLAG  (OUTPUT)  TERMINATION  FLAG. . .POSSIBLE  VALUES  ARE  GR018130 

0 NORMAL  COMPLETION,  E SATISFIES  GR018140 


c 

E<EPS  AND  E<EPS*DABS(R) 

GR018150 

c 

1 

NORMAL  COMPLETION,  E SATISFIES 

GR018160 

c 

E<EPS,  BUT  NOT  RELATIVE  ERROR  REQUEST 

GR018170 

c 

2 

NORMAL  COMPLETION,  E SATISFIES 

GR018180 

c 

E<EPS*DA8S(R) 

GR018190 

c 

3 

NORMAL  COMPLETION  BUT  EPS  WAS  TOO  SMALL 

TO 

GR018200 

c 

SATISFY  ABSOLUTE  OR  RELATIVE  ERROR  REQUEST. 

GR018210 

c 

4 

ABORTED  CALCULATION  BECAUSE  OF  SERIOUS 

ROUNDING 

GR018220 

c 

ERROR.  PROBABLY  E AND  R ARE  CONSISTENT. 

GR018230 

c 

5 

ABORTED  CALCULATION  BECAUSE  OF  INSUFFICIENT  STORAGR018240 

c 

R AND  E ARE  CONSISTENT.  PERHAPS  INCREASING  NMAGR018250 

c 

WILL  PRODUCE  BETTER  RESULTS. 

GRO 18260 

c 

6 

ABORTED  CALCULATION  BECAUSE  OF  SERIOUS 

DIFFICULTIGR018270 

c 

MEETING  YOUR  ERROR  REQUEST. 

GR018280 

c 

7 

ABORTED  CALCULATION  BECAUSE  EITHER  EPS, 

NINT  OR 

NGR018290 

c 

HAS  BEEN  SET  TO  AN  ILLEGAL  VALUE. 

GR018300 

c 

8 

ABORTED  CALCULATION  BECAUSE  YOU  SET  NINT>1  BUT 

FOGR018310 

c 

TO  SET  W( 1 , 1 )=A  AND  W(NINT+1 , 1 )=B 

GR018320 

c 

GR018330 

c 

TYPICAL 

PROBLEM  SET  UP 

GR018340 

c 

GR018350 

c 

DIMENSION  U(50,6) 

GR018360 

c 

LOGICAL  RST 

GRO 18370 

c 

EXTERNAL  F 

GR018380 

c 

A=0.0 

GR018390 

c 

B=1 .0 

GR018400 

c 

W( 1 , 1 )=A 

GR018410 

c 

W(2,1 )=.3 

[SET  INTERNAL  PARTITION  POINT  AT  .3] 

GR018420 

c 

W(3, 1 )=B 

GR018430 

c 

N I NT =2 

[INITIAL  PARTITION  HAS  2 INTERVALS] 

GR018440 

c 

RST=. FALSE. 

GR018450 

c 

EPS=.001 

GR018460 

c 

NMAX=50 

GR018470 

c 

GR018480 

c 

CALL  DGLQ1(F,A, 

B,EPS,R, E,NINT,RST,W,NMAX, FMIN, FMAX,KF, I FLAG) 

GR018490 

c 

GR018500 

c 

1 A,8,EPS 

,R,E,NINT,FMIN, FMAX,KF, I FLAG 

GR018510 

c 

IF(EPS.EQ.  .0001  .OR.  I FLAG.GE.3)STOP 

GR018520 

c 

RST=.TRUE 

GR018530 

c 

EPS=.0001 

[ASK  FOR  ANOTHER  DIGIT] 

GR018540 

c 

GO  TO  1 

GR018550 

c 

END 

GRO 18560 

c 

FUNCTION  F(X) 

GR018570 

c 

IF(X.LT.  .3) 

GR018580 

c 

1 THEN 

GR018590 

c 

F=X**(0.2D0)*ALOG(X) 

GR018600 

c 

ELSE 

GR018610 

c 

F=DSIN(X) 

GR018620 

c 

END  IF 

GR018630 

c 

RETURN 

GR018640 

c 

END 

GR018650 

c 

GR018660 

c 

END  OF  DO 

CUMENTATION 

GRO 18670 

c 

GR018680 

INTEGER  C 

GR018690 

DOUBLE  PRECISION  A , B , E , EB , EPMACH , EPS , FMAX , FMAXL , FMAXR , FM I N , FM I NL 

1 , FMINR, FMN, FMX,Rf RAB, RABS,RAV,SIGN, T,TE 

2 ,TE1 , TE2,TR, TR1 , TR2, UFLOW, W,XM, F 
DIMENSION  W(NMAX,6) 

EXTERNAL  F 

GR018700 
GR018710 
GRO 18720 
GR018730 
GR018740 

LOGICAL  RST.DEBUG 

GR018750 

DATA  SIGN  /-I. DO/ 

GRO 18760 

C THE  FOLLOWING  DATA  ARE  FOR  PERKIN-ELMER  COMPUTERS 

GRO 18770 

DATA  EPMACH , UFL0W/Z341 0000000000000 , Z001 000000000000/ 
C THE  FOLLOWING  DATA  ARE  FOR  VAX  COMPUTERS 

GR018780 
GRO 18790 

C DATA  EPMACH  f UFLOW/Z0000250000000000 , Z0000008000000000/ 

C 

GR018800 

GR018810 

C EPMACH  = Y1 341 0000000000000' 

GR018820 

C UFLOW  = Y'00100000' 

GR018830 

C EPMACH  = DBLE(16**(- 13)) 

C UFLOW  = DBLE(16**(-65)) 

I F(A.EQ.B)  THEN 
R=0. 

GRO 18840 
GR018850 
GR018860 
GRO 18870 

E=0. 

GR018880 

NINT=0 

GRO 18890 

I FLAG=0 

GR018900 

KF=1 

GR018910 

FMIN=F(A) 

GR018920 

FMAX=FMIN 

GRO 18930 

GO  TO  20 

GR018940 

END  IF 

GR018950 

IF(RST)  THEN 

IFdFLAG.LT. 3)  GO  TO  15 
GO  TO  20 

GRO 18960 
GRO 18970 
GR018980 

END  IF 

GR018990 

KF=0 

GR019000 

IFCEPS  .LE.  0.  .OR.  NMAX  .LE.  1 .OR.  NINT  .LE.  0)  THEN 
I FLAG=7 

GR019010 

GR019020 

GO  TO  20 

GR019030 

END  IF 

GR019040 

IF(NINT.EQ.I)  THEN 
C 1 THEN 

GR019050 

GR019060 

W( 1 , 1 )=A 
W(2,2)=B 
W(1 ,5)=A 
W(1 ,6)=B 
W(2,5)=A 
W<2,6)=8 

W(1 ,2)=A+(B-A)/2.D0 
W(2,1)=A+(B-A)/2.D0 
NINT=2 

GR019070 
GR019080 
GR019090 
GR019100 
GR0191 10 
GR019120 
GR019130 
GR019140 
GR019150 

ELSE 

GR019160 

I F(W( 1 , 1 ) .NE .A  .AND.  WCNINT+1 , 1 ) .NE.B)  THEN 
IFLAG=8 

GR019170 

GR019180 

GO  TO  20 

GR019190 

END  IF 

GR019200 

W(1 ,5)=A 
DO  89  1=1 , NINT 
W< I ,2)=W< 1+1 , 1 ) 
W( I ,5)=W( 1,1) 

GR019210 

GR019220 

GR019230 

GR019240 

W( I ,6)=W(I ,2) 
89  CONTINUE 

GR019250 

GR019260 

END  IF 

GR 01 9270 

C 

GR019280 

DEBUG=. FALSE. 

GR019290 

I FLAG  = 0 

GR019300 

IROFF=0 

GR019310 

RABS=0.0 

GR019320 

DO  3 1=1 ,NINT 

GR019330 

CALL  GL15T(F,W(If1),W(If2)fW(I,5),W(I,6), 

1 W(I ,4),W(I ,3),RAB,RAV, FMN, FMX) 

I F(DEBUG)  WRITE (7  ,*) 1 INITIALIZE (W( I , J) ,J=1 ,6) 
KF=KF+15 

GRO 19340 
GR019350 
GR019360 
GRO 19370 

IF(I.EQ.I) 
1 THEN 

GR019380 
GRO 19390 

R=W( 1,4) 

E=W( I ,3) 
RABS=RABS+RAB 

GR019400 

GR019410 

GR019420 

FMIN=FMN 

GR019430 

FMAX=FMX 

GR019440 

ELSE 

GR019450 

R=R+W( I ,4) 
E=E+U( I ,3) 
RABS=RABS+RAB 

GR019460 
GR019470 
GRO 19480 

FMAX=DMAX1 ( FMAX , FMX ) 
FMIN=DMIN1 ( FMIN, FMN) 

END  IF 

GR019490 

GR019500 

GR019510 

3 CONTINUE 

GR019520 

DO  10  I=NINT+1,NMAX 

GR019530 

W( I ,3)  = 0. 
10  CONTINUE 

GR019540 

GR019550 

15  CONTINUE 

GR019560 

C 

GR019570 

C MAIN  SUBPROGRAM  LOOP 

GR019580 

C 

GR019590 

I F ( 1 00 . DO*EPMACH*RABS . GE .DABS( R ) .AND.  E.LT.EPSJGO  TO  20 
EB=DMAX1(100.D0*UFLOW,DMAXl (EPS,50.D0*EPMACH)*DABS(R) ) 
IF(E.LE.EB)  GO  TO  20 

GR019600 
GR019610 
GRO 19620 

IF  (NINT.LT.NMAX) 
1 THEN 

GR019630 

GR019640 

NINT  = NINT+1 

GR019650 

C = NINT 

GR019660 

ELSE 

GRO 19670 

C=0 

GR019680 

16  I F(C.EQ.NMAX)  THEN 
I FLAG=5 

GRO 19690 
GR019700 

GO  TO  20 

GR019710 

END  IF 

GRO 19720 

C=C+1 

GR019730 

I F(W(C,3) .GT.O.DO)  GO  TO  16 
END  IF 

GR019740 

GR019750 

C LOC=ISAMAX(NINT,W(1 , 3) , 1 ) 

CALL  LSE(W,NMAX,NINT,LOC) 

I F(DEBUG)WRITE(7,200)LOC, (W(LOC, I),I=1,6),R,E 
XM  = W(LOC, 1 )+(W(LOC,2)-U(LOC, 1 ))/2. 

GR019760 
GR019770 
GR019780 
GRO 19790 

IF  ( (DMAX1 (DA3S(U(L0C , 1 ) ) ,DABS(W(LOC, 2) ) ) ) .GT. 

GR019800 

1 ( ( 1 .D0+100.D0*EPMACH)*(DABS(XM)+0. 1D+04*UFLOW))) 

GR019810 

2 THEN 

GR019820 

CALL  GL15T ( F,W(LOC, 1 ) ,XM, W(L0C,5) ,W(L0C,6) , 

GR019830 

1 TR1 ,TE1 , RAB,RAV, FMINL, FMAXL) 

GR019840 

KF=KF+15 

GR019850 

IF  (TE1.LT. (EB*(XM-U(LOC, 1))/(B-A))) 

GR019860 

A TE1  = TE1*SIGN 

GRO 19870 

CALL  GL15T(F,XMlW(LOC,2),W(LCCt5),U(LOC,6), 

GR019880 

1 TR2,TE2,RA8,RAV, FMINR, FMAXR) 

GR019890 

KF=KF+15 

GR019900 

FMIN=DMIN1(FMIN, FMINL, FMINR) 

GR019910 

FMAX=DMAX1 (FMAX, FMAXL, FMAXR) 

GR019920 

IF  (TE2.LT. (EB*(W(L0C,2)-XM)/(B-A)))  TE2=TE2*SIGN 

GR019930 

TE  = DABS(W(L0C,3)) 

GR019940 

TR  = U(LOC, 4) 

GR019950 

U(C, 3)  = TE2 

GR019960 

W(C,4)  = TR2 

GR019970 

W(C,1)  = XM 

GR019980 

W(C,2)  = W(L0C,2) 

GR019990 

W(C,5)  = W(L0C,5) 

GR020000 

W(C,6)  = W(LOC, 6) 

GR020010 

W(L0C,3)  = TE1 

GR020020 

W(L0C,4)  = TR1 

GR020030 

W(L0C,2)  = XM 

GR020040 

I F(DEBUG)WRITE(7, 200)C, (U(C,K),K=1,6) 

GR020050 

I F( DEBUG) WRITE (7, 200) LOC, (W(L0C,K) ,K=1 ,6) 

GR020060 

E = E-TE+(DABS(TE1)+DABS(TE2)) 

GR020070 

R = R-TR+(TR1+TR2) 

GR020080 

I F (DEBUG)WRITE(7  ,*)NINT,R,E 

GR020090 

I F(DABS(DABS(TE1 )+DABS(TE2)-TE) .LT. .001D0*TE)  THEN 

GR020100 

IROFF=IROFF+1 

GR0201 10 

I F( IROFF.GE .10)  THEN 

GR020120 

I FLAG=4 

GR020130 

GO  TO  20 

GR020140 

END  IF 

GR020150 

END  IF 

GR020160 

ELSE 

GR020170 

IF  (EB.GT.U(LOC,3)) 

GR020180 

1 THEN 

GR020190 

W(LOC,3)  = 0. 

GR020200 

ELSE 

GR020210 

IFLAG=6 

GR020220 

GO  TO  20 

GR020230 

END  IF 

GR020240 

END  IF 

GR020250 

GO  TO  15 

GR020260 

GR020270 

ALL  EXITS  FROM  HERE 

GR020280 

GR020290 

20  CONTINUE 

GR020300 

I F ( I FLAG . GE . 4 )RETURN 

GR020310 

I FLAG=3 

GR020320 

T=EPS*DABS(R) 

GR020330 

I F(E.GT.EPS  .AND.  E.GT.T)RETURN 

GR020340 

I FLAG=2 

I F(E.GT.EPS  .AND.  E.LT.T)RETURN 
IFLAG=1 

IF(E.LT.EPS  .AND.  E.GT.T)RETURN 

IFLAG=0 

RETURN 

200  FORMAT  ( I4,8(E1 1 .3)) 

END 


SUBROUTINE  LSE  (WORK,NMAX,NI , LOC) 

C 

C 

THIS  SUBPROGRAM  FINDS  THE  CELL  IN  THE  WORK  AREA 
THAT  HAS  THE  LARGEST  ABSOLUTE  ERROR  AND 
RETURNS  THE  LOCATION  OF  THAT  CELL. 


DOUBLE  PRECISION  ERROR, WORK 
DIMENSION  WORK(NMAX,6) 

INITIALIZE  VARIABLES 

ERROR  = DABS(WORK( 1 ,3) ) 

LOC  = 1 

MAIN  SUBPROGRAM  LOOP 
DO  20  1=1, NI 

IF  (DABS(WORK(I,3)).GT.ERROR)  THEN 
ERROR  = DABS(WORK( 1,3)) 

LOC  = I 
END  IF 
20  CONTINUE 
RETURN 
END 


SUBROUTINE  GL15T(F,A,B,XL,XR,R,AE,RA, 
1 RASC, FMIN, FMAX) 


C***AUTHORS 

C 

C 


ROBERT  PIESSENS  AND  ELISE  DE  DONCKER 
APPL.  MATH.  AND  PROGR.  DIV.-  K.U. LEUVEN 
DAVID  KAHANER,  NBS  WASHINGTON 


PURPOSE 

TO  COMPUTE  I = INTEGRAL  OF  G(X)  OVER  (A, B) , 
WITH  ERROR  ESTIMATE 
J = INTEGRAL  OF  DABS(G)  OVER  (A,B) 


PARAMETERS 
ON  ENTRY 

F - DOUBLE  PRECISION 


GR020350 

GR020360 

GR020370 

GR020380 

GR020390 

GR020400 

GR020410 

GR020420 

GR020430 

GR020440 

GR020450 

GR020460 

GR020470 

GR020480 

GR020490 

GR020500 

GR020510 

GR020520 

GR020530 

GR020540 

GR020550 

GR020560 

GR020570 

GR020580 

GR020590 

GR020600 

GR020610 

GR020620 

GR020630 

GR020640 

GR020650 

GR020660 

GR020670 

GR020680 

GR020690 

GR020700 

GR020710 

GR020720 

GR020730 

GR020740 

GR020750 

GR020760 

GR020770 

GR020780 

GR020790 

GR020800 

GR020810 

GR020820 

GR020830 

GR020840 

GR020850 

GR020860 

GR020870 

GR020880 

GR020890 


c 

FUNCTION  SUBPROGRAM  DEFINING  THE  INTEGRAND 

GR020900 

c 

FUNCTION  F(X) . THE  ACTUAL  NAME  FOR  F NEEDS 

GR02C910 

c 

TO  BE  DECLARED  E X T E R N A L IN  THE 

GR020920 

c 

CALLING  PROGRAM. 

GR020930 

c 

THE  FUNCTION  G(X)  IS  DEFINED  TO  BE 

GR020940 

c 

G(X)=F(PHI(X))*PHIP(X) 

GR020950 

c 

WHERE  PHI (X)  IS  THE  CUBIC  GIVEN  BY 

GR020960 

c 

THE  ARITHMETIC  STATEMENT  FUNCTION  BELOW. 

GR020970 

c 

PHIP(X)  IS  ITS  DERIVATIVE.  THE  VARIABLES 

GR020980 

c 

XL  AND  XR  ARE  THE  LEFT  AND  RIGHT  ENDPOINTS 

GR020990 

c 

OF  A PARENT  INTERVAL  OF  WHICH  (A,B)  IS  A PART. 

GR021000 

c 

GR021010 

c 

A - DOUBLE  PRECISION 

GR021020 

c 

LOWER  LIMIT  OF  INTEGRATION 

GR021030 

c 

GR021040 

c 

B - DOUBLE  PRECISION 

GR021050 

c 

UPINTG  LIMIT  OF  INTEGRATION 

GR021060 

c 

GR021070 

c 

XL  - DOUBLE  PRECISION 

GR021080 

c 

LEFT  ENDPOINT  OF  PARENT  INTERVAL 

GR021090 

c 

GR021100 

c 

XR  - DOUBLE  PRECISION 

GR021110 

c 

RIGHT  ENDPOINT  OF  PARENT  INTERVAL 

GR021120 

c 

GR021130 

c 

ON  RETURN 

GR021140 

c 

R - DOUBLE  PRECISION 

GR021150 

c 

APPROXIMATION  TO  THE  INTEGRAL  I 

GR021160 

c 

R IS  COMPUTED  BY  APPLYING  THE  15-POINT 

GR021170 

c 

KRONROD  RULE  (RESK)  OBTAINED  3Y  OPTIMAL 

GR021 180 

c 

ADDITION  OF  ABSCISSAE  TO  THE  7-POINT  GAUSS 

GR021190 

c 

RULE  (RESG) . 

GR021200 

c 

GR021210 

c 

AE  - DOUBLE  PRECISION 

GR021220 

c 

ESTIMATE  OF  THE  MODULUS  OF  THE  ABSOLUTE  ERROR, 

GR021230 

c 

WHICH  SHOULD  NOT  EXCEED  DABS(I-R) 

GR021240 

c 

GR021250 

c 

RA  - DOUBLE  PRECISION 

GR021260 

c 

APPROXIMATION  TO  THE  INTEGRAL  J 

GR021270 

c 

GR021280 

c 

RASC  - DOUBLE  PRECISION 

GR021290 

c 

APPROXIMATION  TO  THE  INTEGRAL  OF  DA3S(G- I/(B-A) )GR021300 

c 

OVER  (A, B) 

GR021310 

c 

GR021320 

c 

FMAX,  FMIN  - DOUBLE  PRECISION 

GR021330 

c 

MAX  AND  MIN  VALUES  OF  THE  FUNCTION  F ON  (A,B) 

GR021340 

c 

SUBROUTINES  OR  FUNCTIONS  NEEDED 

GR021350 

c 

- F (USER-PROVIDED  FUNCTION) 

GR021360 

c 

- R1MACH 

GR021370 

c 

- FORTRAN  DABS,  DMAX1 , DMIN1 

GR021380 

c 

GR021390 

c 

. GR021400 

C***END 

PROLOGUE 

GR021410 

c 

GR021420 

DOUBLE  PRECISION  A,AE,B,DHLGTH,EPMACH, F, FC, FSUM, FVAL1 , FVAL2,  GR021430 

* FV1 , FV2,HLGTH,RA,RASC,RESG,RESK,RESKH,R,R1MACH,UFLCW,  GR021440 


* WG,WGK,XGK 

DOUBLE  PRECISION  XL, XR, CENTR, ABSC,U, FMAX, FMIN,PHI , PHIP 
INTEGER  J , JTW, JTWM1 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


C 


DIMENSION  FV1(7),FV2(7),WG(4),WGK(8),XGK(8) 


THE  ABSCISSAE  AND  WEIGHTS  ARE  GIVEN  FOR  THE  INTERVAL  (-1,1) 
BECAUSE  OF  SYMMETRY  ONLY  THE  POSITIVE  ABSCISSAE  AND  THEIR 
CORRESPONDING  WEIGHTS  ARE  GIVEN. 


XGK  - ABSCISSAE  OF  THE  15-POINT  KRONROD  RULE 

XGK(2) , XGK(4) , ...  ABSCISSAE  OF  THE  7-POINT 
GAUSS  RULE 

XGK<1),  XGK(3) , ...  ABSCISSAE  WHICH  ARE  OPTIMALLY 
ADDED  TO  THE  7-POINT  GAUSS  RULE 

WGK  - WEIGHTS  OF  THE  15-POINT  KRONROD  RULE 

WG  - WEIGHTS  OF  THE  7-POINT  GAUSS  RULE 


DATA  XGK( 1 ) ,XGK(2) ,XGK(3) ,XGK(4) ,XGK(5) ,XGK(6) ,XGK(7) ,XGK(8)/ 
0 .991455371 1 2081 26D+00 , 0 . 9491 0791 23427585D+00 , 

0 . 8648644233597691 D+00 , 0.7415311 855993944D+00 , 

0 . 586087235467691 1 D+00 , 0 . 405845 1 5 1 3773972D+00 , 

0.2077849550078985D+00,  0.0D+00  / 

DATA  WGK( 1 ) ,WGK(2) ,WGK(3) ,WGK(4) ,WGK(5) , WGK(6) ,WGK(7) ,WGK(8)/ 
0 . 2293532201 052922D - 01 , 0 . 6309209262997855D - 01 , 

0. 1047900103222502D+00,  0. 14065325971552590+00, 

0 . 1 6900472 66392679D+00 , 0.1 903505780647854D+00 , 


* 0.2044329400752989D+00, 
DATA  WG(1),WG(2),WG(3),WG(4)/ 

* 0.1 294849661 688697D+00, 

* 0.3818300505051 1890+00, 


0.20948214108472780+00/ 

0 . 279705391 4892767D+00 , 
0.41 79591 836734694D+00/ 


DATA  EPMACH , UFL0W/Z341 0000000000000 , Z001 0000000000000/ 


PHI (U)=XR- (XR-XL)*U*U*(2.D0*U+3.D0) 
PHIP(U)=-6.D0*U*(U+1 .DO) 

C 


c 

c 

LIST  OF 

MAJOR  VARIABLES 

c 

c 

CENTR 

- MID  POINT  OF  THE  INTERVAL 

c 

HLGTH 

- HALF-LENGTH  OF  THE  INTERVAL 

c 

ABSC 

- ABSCISSA 

c 

FVAL* 

- FUNCTION  VALUE 

c 

RESG 

- R OF  THE  7-POINT  GAUSS  FORMULA 

c 

RESK 

- R OF  THE  15-POINT  KRONROD  FORMULA 

c 

c 

RESKH 

- APPROXIMATION  TO  THE  MEAN  VALUE  OF  F OVER  (A,B) 
I.E.  TO  I/CB-A) 

c 

c 

MACHINE 

DEPENDENT  CONSTANTS 

C 

C EPMACH  IS  THE  LARGEST  RELATIVE  SPACING. 

C UFLOW  IS  THE  SMALLEST  POSITIVE  MAGNITUDE. 


GR021450 

GR021460 

GR021470 

GR021480 

GR021490 

GR021500 

GR021510 

GR021520 

GR021530 

GR021540 

GR021550 

GR021560 

GR021570 

GR021580 

GR021590 

GR021600 

GR021610 

GR021620 

GR021630 

GR021640 

GR021650 

GR021660 

GR021670 

GR021680 

GR021690 

GR021700 

GR021710 

GR021720 

GR021730 

GR021740 

GR021750 

GR021760 

GR021770 

GR021780 

GR021790 

GR021800 

GR021810 

GR021820 

GR021830 

GR021840 

GR021850 

GR021860 

GR021870 

GR021880 

GR021890 

GR021900 

GR021910 

GR021920 

GR021930 

GR021940 

GR021950 

GR021960 

GR021970 

GR021980 

GR021990 


c 

GR022000 

C***F IRST  EXECUTABLE  STATEMENT 

GR022010 

C 

EPMACH  = DBLE( 16**( - 13) ) 

GR022020 

C 

UFLOW  = DBLE(16**( -65)) 

GR022030 

c 

GR022040 

HLGTH  = 0.5D+00*(B-A) 

GR022050 

CENTR  = A+HLGTH 

GR022060 

DHLGTH  = OABS(HLGTH) 

GR022070 

c 

GR022080 

c 

COMPUTE  THE  15-POINT  KRONROD  APPROXIMATION  TO 

GR022090 

c 

THE  INTEGRAL,  AND  ESTIMATE  THE  ABSOLUTE  ERROR. 

GR022100 

c 

GR0221 10 

U= ( CENTR -XR)/(XR- XL) 

GR022120 

FMIN=F(PHI (U) ) 

GR022130 

FMAX=FMI N 

GR022140 

FC=FMIN*PHIP(U) 

GR022150 

RESG  = FC*UG(4) 

GR022160 

RESK  = FC*WGK(8) 

GR022170 

RA  = DABS(RESK) 

GR022180 

DO  10  J=1 ,3 

GR022190 

JTU  = J*2 

GR022200 

ABSC  = HLGTH*XGK( JTU) 

GR022210 

U=( CENTR -ABSC-XR)/(XR- XL) 

GR022220 

FVAL1  = F(PHI(U)) 

GR022230 

FMAX=DMAX 1 ( FMAX , FVAL1 ) 

GR022240 

FMIN=DMIN1 (FMIN, FVAL1 ) 

GR022250 

FVAL1=FVAL1*PHIP(U) 

GR022260 

U=(CENTR+ABSC-XR)/(XR-XL) 

GR022270 

FVAL2=F(PHI (U) ) 

GR022280 

FMAX=DMAX 1 ( FMAX , FVAL2) 

GR022290 

FMIN=DMIN1 (FMIN, FVAL2) 

GR022300 

FVAL2=FVAL2*PHIP(U) 

GR022310 

FVI(JTW)  = FVAL1 

GR022320 

FV2CJTW)  = FVAL2 

GR022330 

FSUM  = FVAL1+FVAL2 

GR022340 

RESG  = RESG+WGC J)*FSUM 

GR022350 

RESK  = RESK+WGK( JTU)* FSUM 

GR022360 

RA  = RA+WGK( JTU)*(DABS( FVAL1 )+DABS( FVAL2) ) 

GR022370 

10  CONTINUE 

GR022380 

DO  15  J = 1,4 

GR022390 

JTWM1  = J*2- 1 

GR022400 

ABSC  = HLGTH*XGK( JTWM1 ) 

GR022410 

U= ( CENTR - ABSC -XR)/ (XR-XL) 

GR022420 

FVAL1=F(PHI <U) ) 

GR022430 

FMAX=DMAX1 ( FMAX, FVAL1 ) 

GR022440 

FMIN=DMIN1 (FMIN, FVAL1 ) 

GR022450 

FVAL1=FVAL1*PHIP(U) 

GR022460 

U=(CENTR+ABSC-XR)/( XR-XL) 

GR022470 

FVAL2=F(PHI(U)) 

GR022480 

FMAX=DMAX1 ( FMAX, FVAL2) 

GR022490 

FMI N=DMIN1 ( FMIN , FVAL2) 

GR022500 

FVAL2=FVAL2*PHIP(U) 

GR022510 

FVI(JTWMI)  = FVAL1 

GR022520 

FV2( JTWM1 ) = FVAL2 

GR022530 

FSUM  = FVAL1+FVAL2 

GR022540 

RESK  = RESK+WGK( JTWM1 )*FSUM 

RA  = RA+WGK( JTWM1 )*(DABS( FVAL1 )+OABS( FVAL2) ) 

15  CONTINUE 

RESKH  = RESK*0.5D+00 

RASC  = WGK(8)*DABS(FC-RESKH) 

DO  20  J=1 ,7 

RASC  = RASC+WGK( J )*(DABS( FV1 ( J ) - RESKH )+DABS( FV2( J ) - RESKH ) ) 
20  CONTINUE 

R = RESK*HLGTH 

RA  = RA*DHLGTH 

RASC  = RASC*DHLGTH 

AE  = DABS( (RESK-RESG)*HLGTH) 

I F ( RASC . NE . 0 . 0D+00 . AND . AE . NE . 0 . 0D+00 ) 

* AE  = RASC’DMINI (0. 1D+01 , 

* (0.2D+03*AE/RASC)**1 .5D+00) 

I F(RA.GT.UFLOW/(0.5D+02*EPMACH) ) AE  = DMAX1 

* ( (EPMACH*0.5D+02)*RA,AE) 

RETURN 

END 


C 

SUBROUTINE  QUADRT ( A , ROOT 1 , ROOT2 , TOL , NOROOT ) 
C 

C 

C TO  SOLVE  ANY  QUADRT I C EQUATION  OF  THE  FORM 
C 

C A(1 )*X**2+A(2)*X+A(3)  = 0 

C 

DOUBLE  PRECISION  A, ROOT  1 , ROOT2,X,Y,Z,U,ZERO 
DIMENSION  A(3) ,ROOT1 (2) , ROOT2L2) 

C 

ZERO  = TOL/IO.O 
C 

IF  (DABS(A(1 )) -ZERO) 2, 2, 1 

2 IF  (DABS(A(2) )-ZERO)  3,3,4 

3 NOROOT  =0 
RETURN 

C 

1 NOROOT = 2 

X = A(2)*A(2)  - 4.0*A( 1 )*A(3) 

Y = A( 1 ) + A( 1 ) 

Z = DSQRT (DA8S(X) )/Y 
U = -A(2)/Y 
C 

IF  (X.LT.0.0)  GO  TO  7 
ROOT1 ( 1 )=W+Z 
ROOT1 (2)=0.0 
ROOT2( 1 )=W-Z 
ROOT2(2)=0.0 
RETURN 
C 

7 ROOT1 ( 1 )=W 
ROOT1 (2)=Z 
ROOT2( 1 )=W 


GR022550 

GR022560 

GR022570 

GR022580 

GR022590 

GR022600 

GR022610 

GR022620 

GR022630 

GR022640 

GR022650 

GR022660 

GR022670 

GR022680 

GR022690 

GR022700 

GR022710 

GR022720 

GR022730 

GR022740 

GR022750 

GR022760 

GR022770 

GR022780 

GR022790 

GR022S00 

GR022810 

GR022820 

GR022830 

GR022840 

GR022850 

GR022860 

GR022870 

GR022880 

GR022890 

GR022900 

GR022910 

GR022920 

GR022930 

GR022940 

GR022950 

GR022960 

GR022970 

GR022980 

GR022990 

GR023000 

GR023010 

GR023020 

GR023030 

GR023040 

GR023050 

GR023060 

GR023070 

GR023080 

GR023090 


ROOT2(2)=-Z 

RETURN 

C 

4 N0R00T=1 

ROOTI(I)  = -A(3)/A(2) 

ROOT1 (2)  = 0.0 

RETURN 

END 

£•»***************•*************************************************** 


c 

c** 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE  DCRTNI (X, F,DERF,XST,EPS, IEND, IER) 

★ **★*★★*★★*★★★*★★★*★★*★★★★★★★★★★★**★★**★★*★*★**★**★***★★*★**★**** 

PURPOSE 

TO  SOLVE  PHI  = 0 

BY  MEANS  OF  NEUTON-S  ITERATION  METHOD. 

USAGE 

CALL  DCRTNI  (X, F, DERF, XST, EPS, IEND, IER) 

DESCRIPTION  OF  PARAMETERS 

X - DOUBLE  COMPLEX  RESULTANT  ROOT  OF  EQUATION  F(X)=0. 

F - DOUBLE  COMPLEX  RESULTANT  FUNCTION  VALUE  AT 

ROOT  X. 

DERF  - DOUBLE  COMPLEX  RESULTANT  VALUE  OF  DERIVATIVE 
AT  ROOT  X. 

XST  - DOUBLE  COMPLEX  INPUT  VALUE  WHICH  SPECIFIES  THE 
INITIAL  GUESS  OF  THE  ROOT  X. 

EPS  - SINGLE  PRECISION  INPUT  VALUE  WHICH  SPECIFIES  THE 
UPPER  BOUND  OF  THE  ERROR  OF  RESULT  X. 

IEND  - MAXIMUM  NUMBER  OF  ITERATION  STEPS  SPECIFIED. 

IER  - RESULTANT  ERROR  PARAMETER  CODED  AS  FOLLOWS 
IER=0  - NO  ERROR, 

IER=1  • NO  CONVERGENCE  AFTER  IEND  ITERATION  STEPS 
IER=2  - AT  ANY  ITERATION  STEP  DERIVATIVE  DERF  WAS 
EQUAL  TO  ZERO. 


SUBROUTINES  AND  FUNCTION  SUBPROGRAMS  REQUIRED 
PHIO.AND  DPHIO 

METHOD 

SOLUTION  OF  EQUATION  F(X)=0  IS  DONE  BY  MEANS  OF  NEWTON-S 
ITERATION  METHOO,  WHICH  STARTS  AT  THE  INITIAL  GUESS  XST  OF 
A ROOT  X.  CONVERGENCE  IS  QUADRATIC  IF  THE  DERIVATIVE  OF 
F(X)  AT  ROOT  X IS  NOT  EQUAL  TO  ZERO.  ONE  ITERATION  STEP 
REQUIRES  ONE  EVALUATION  OF  F(X)  AND  ONE  EVALUATION  OF  THE 
DERIVATIVE  OF  F<X). 

FOR  REFERENCE,  SEE  R.  ZURMUEHL,  PRAKTISCHE  MATHEMATIK  FUER 
INGENIEURE  UND  PHYSIKER,  SPRINGER,  BERLIN/GCETTINGEN/ 
HEIDELBERG,  1963,  PP.12-1 7. 


DOUBLE  COMPLEX  X,F, DERF, XST, TOL,DX,PHIO, DPHIO 


GR023100 
GR0231 10 
GR023120 
GR023130 
GR023140 
GR023150 
GR023160 
GR023170 
GR023180 
GR023190 
GR023200 
GR023210 
GR023220 
GR023230 
GR023240 
GR023250 
GR023260 
GR023270 
GR023280 
GR023290 
GR023300 
GR023310 
GR023320 
GR023330 
GR023340 
GR023350 
GR023360 
GR023370 
GR023380 
GR023390 
GR023400 
GR023410 
GR023420 
GR023430 
, GR023440 
GR023450 
GR 023460 
GR023470 
GR023480 
GR023490 
GR023500 
GR023510 
GR023520 
GR023530 
GR023540 
GR023550 
GR023560 
GR023570 
GR023580 
GR023590 
GR023600 
GR 0236 10 
GR023620 
GR023630 
GR023640 


DOUBLE  PRECISION  A,DTOL,DTOLF 
EXTERNAL  PHI0,DPHI0 
C PREPARE  ITERATION 
IER=0 
X=XST 
TOL=X 

F = PHIO(TOL) 

DERF  = DPHIO(TOL) 

DTOLF=100.D0*EPS 

C 

C 

C START  ITERATION  LOOP 

DO  6 1=1 , IEND 
I F(CDABS(F))1 ,7, 1 
C 

C EQUATION  IS  NOT  SATISFIED  BY  X 

1 I F(CDABS(DERF))2,8,2 
C 

C ITERATION  IS  POSSIBLE 

2 DX=F/DERF 
X=X-DX 
TOL=X 

F = PHIO(TOL) 

DERF  = DPHIO(TOL) 

C 

C TEST  ON  SATISFACTORY  ACCURACY 

DTOL=DBLE(EPS) 

A=CDABS(X) 

I F(A- 1 .D0)4, 4,3 

3 DTOL=DTOL*A 

4 I F(CDABS(DX)-DTOL)5,5,6 

5 IF(CDABS(F)-DT0LF)7,7,6 

6 CONTINUE 

C END  OF  ITERATION  LOOP 

C 

C 

C NO  CONVERGENCE  AFTER  IEND  ITERATION  STEPS.  ERROR  RETURN. 
IER=1 

7 RETURN 
C 
C 

8 


SUBROUTINE  RSORT(A,B,R,N,M,MS) 

C 

DIMENSION  A(1),B(1),R(1) 

C 

C MOVE  SORTING  KEY  VECTOR  TO  FIRST  COLUMN  OF  OUTPUT  MATRIX 

C AND  3UILD  ORIGINAL  SEQUENCE  LIST  IN  SECOND  COLUMN 

C 

DO  10  1=1, N 


ERROR  RETURN  IN  CASE  OF  ZERO  DIVISOR 

IER=2 

RETURN 

END 


GR023650 

GR023660 

GR023670 

GR023680 

GR023690 

GR023700 

GR023710 

GR023720 

GR023730 

GR023740 

GR023750 

GR023760 

GR023770 

GR023780 

GR023790 

GR023800 

GR023810 

GR023820 

GR023830 

GR023840 

GR023850 

GR023860 

GR023870 

GR023880 

GR023890 

GR023900 

GR023910 

GR023920 

GR023930 

GR023940 

GR023950 

GR023960 

GR023970 

GR023980 

GR023990 

GR024000 

GR024010 

GR024020 

GR024030 

GR024040 

GR024050 

GR024060 

GR024070 

GR024080 

GR024090 

GR024100 

GR024110 

GR024120 

GR024130 

GR024140 

GR024150 

GR024160 

GR 024 170 

GR024180 

GR024190 


R(  I )=B( I ) 

GR024200 

I2=I>M 

GR024210 

10 

R(I2)=I 

GR024220 

c 

GR024230 

c 

SORT  ELEMENTS  IN  SORTING  KEY  VECTOR 

(ORIGINAL  SEQUENCE  LIST 

GR024240 

r 

IS  RESEQUENCED  ACCORDINGLY) 

GR024250 

C 

GR024260 

L=N+1 

GR024270 

20 

I S0RT=0 

GR024280 

L=L- 1 

GR024290 

I F(L  .LT.  2)  GO  TO  50 

GR024300 

DO  40  1=2, L 

GR024310 

I F(R( I)-R(I-I))  30,40,40 

GR024320 

30 

ISORT=1 

GR024330 

RSAVE=R< I ) 

GR024340 

R( I )=R( 1*1) 

GR024350 

R( I - 1 )=RSAVE 

GR024360 

I2=I+N 

GR024370 

SAVER=R( 12) 

GR024380 

R< 1 2)=R( 12-1 ) 

GR024390 

R ( 1 2 - 1 )=SAVER 

GR024400 

40 

CONTINUE 

GR024410 

I F( I SORT)  20,50,20 

GR024420 

c 

GR024430 

c 

MOVE  ROWS  FROM  MATRIX  A TO  MATRIX 

R 

(NUMBER  IN  SECOND  COLUMN 

GR024440 

c 

OF  R REPRESENTS  ROW  NUMBER  OF  MATRIX 

A TO  BE  MOVED) 

GR024450 

c 

GR024460 

50 

DO  80  1=1, N 

GR024470 

c 

GR024480 

c 

GET  ROW  NUMBER  IN  MATRIX  A 

GR024490 

c 

GR024500 

I2=I+N 

GR024510 

I N=R ( I 2 ) 

GR024520 

c 

GR024530 

IR=I -N 

GR024540 

DO  80  J=1,M 

GR024550 

c 

GR024560 

c 

LOCATE  ELEMENT  IN  OUTPUT  MATRIX 

GR024570 

c 

GR024580 

IR=IR+N 

GR024590 

c 

GR024600 

c 

LOCATE  ELEMENT  IN  INPUT  MATRIX 

GR024610 

c 

GR024620 

CALL  LOC( IN , J , IA, N,M,MS) 

GR024630 

c 

GR024640 

c 

TEST  FOR  ZERO  ELEMENT  IN  DIAGONAL 

MATRIX 

GR024650 

c 

GR024660 

I FC IA)  60,70,60 

GR024670 

c 

GR024680 

c 

MOVE  ELEMENT  TO  OUTPUT  MATRIX 

GR024690 

c 

GR024700 

60 

R( IR)=A( IA) 

GR024710 

GO  TO  80 

GR024720 

70 

R( I R ) =0 

GR024730 

80 

CONTINUE 

GR024740 

RETURN 

END 

Q- ******************************************************************* 

C 

SUBROUTINE  LOC( I , J, IR,N,M,MS) 

C 


PURPOSE 

COMPUTE  A VECTOR  SUBSCRIPT  FOR  AN  ELEMENT  IN  A MATRIX  OF 
SPECIFIED  STORAGE  MODE 

USAGE 

CALL  LOC  (I,J, IR,N,M,MS) 

DESCRIPTION  OF  PARAMETERS 

I - ROW  NUMBER  OF  ELEMENT 

J - COLUMN  NUMBER  OF  ELEMENT 

IR  - RESULTANT  VECTOR  SUBSCRIPT 

N - NUMBER  OF  ROWS  IN  MATRIX 

M * NUMBER  OF  COLUMNS  IN  MATRIX 

MS  - ONE  DIGIT  NUMBER  FOR  STORAGE  MODE  OF  MATRIX 

0 - GENERAL 

1 - SYMMETRIC 

2 - DIAGONAL 

REMARKS 

NONE 

M IS  UNUSED  BUT  IS  REQUIRED  TO  BE  PRESENT 

SUBROUTINES  AND  FUNCTION  SUBPROGRAMS  REQUIRED 
NONE 


METHOD 

MS=0 

MS=1 


MS=2 


SUBSCRIPT  IS  COMPUTED  FOR  A MATRIX  WITH  N*M  ELEMENTS 
IN  STORAGE  (GENERAL  MATRIX) 

SUBSCRIPT  IS  COMPUTED  FOR  A MATRIX  WITH  N*(N+1)/2  IN 
STORAGE  (UPPER  TRIANGLE  OF  SYMMETRIC  MATRIX).  IF 
ELEMENT  IS  IN  LOWER  TRIANGULAR  PORTION,  SUBSCRIPT  IS 
CORRESPONDING  ELEMENT  IN  UPPER  TRIANGLE. 

SUBSCRIPT  IS  COMPUTED  FOR  A MATRIX  WITH  N ELEMENTS 
IN  STORAGE  (DIAGONAL  ELEMENTS  OF  DIAGONAL  MATRIX). 

IF  ELEMENT  IS  NOT  ON  DIAGONAL  (AND  THEREFORE  NOT  IN 
STORAGE),  IR  IS  SET  TO  ZERO. 


IX=I 

JX=J 

IF(MS-I)  10,20,30 
10  IRX=N*( JX- 1 )+IX 
GO  TO  36 

20  IF(IX-JX)  22,24,24 
22  IRX=IX+( JX*JX- JX)/2 
GO  TO  36 


GR024750 
GR024760 
GR024770 
GR024780 
GR024790 
GR024800 
GR024810 
GR024820 
GR024830 
GR024840 
GR024850 
GR024860 
GR024870 
GR024880 
GR024890 
GR024900 
GR024910 
GR024920 
GR024930 
GR024940 
GR024950 
GR024960 
GR024970 
GR024980 
GR024990 
GR025000 
GR025010 
GR025020 
GR025030 
GR025040 
GR025050 
GR025060 
GR025070 
GR025080 
GR025090 
GR025100 
GR0251 10 
GR025120 
GR025130 
GR025140 
GR025150 
GR025160 
GR025170 
GR025180 
GR025190 
, GR025200 
GR025210 
GR025220 
GR025230 
GR025240 
GR025250 
GR025260 
GR025270 
GR025280 
GR025290 


24 

IRX=JX+( I X* I X - 1 X )/2 

GR025300 

GO  TO  36 

GR025310 

30 

IRX=0 

GR025320 

IF(IX-JX)  36,32,36 

GR025330 

32 

IRX=IX 

GR025340 

36 

IR=IRX 

GR025350 

RETURN 

GR025360 

END 

GR025370 
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